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Introduction 



A deep understanding of the interaction between matter and radiation (including 
electrons and light) is a key issue in order to describe the physical nature and the 
properties of materials. This can be achieved with a joint effort of numerical sim- 
ulations and experiments. In fact, the physical origin of the experimental spectral 
features can often be understood unambiguously with the help of numerical simula- 
tions. 

Nowadays numerical computation of ground state properties of condensed mat- 
ter systems can be successfully treated within the density functional theory (DFT). 
In this context, the problem of solving the Schrodinger equation for the ground state 
of a many body system can be exactly recast into the variational problem of mini- 
mizing an energy functional with respect to the charge density. The success of this 
approach has been shown during the last years by ab initio calculations describing 
the ground state properties of realistic systems, in particular in case of reconstructed 
surfaces. However, ground state properties are not enough to describe those exper- 
iments involving excitations of the electronic system. In most cases, an external 
probe modifies the charge distribution of the sample, producing excited states and 
a dynamical rearrangement of the density. 

In the last years, excited state theories providing an overcome of limits of DFT 
have been proposed. A particularly fruitful attempt to go beyond the ground state 
theory is offered by time dependent density functional theory (TDDFT), providing 
an exact reformulation of quantum mechanics in terms of time evolving density. 
Within this theory, the complexity of the problem is confined to the exchange- 
correlation potential V xc , whose analytic form is unknown but several approxima- 
tions are available in literature, giving correct predictions in many realistic cases. 
However, in many interesting applications, the simpliest approximation (as indepen- 
dent particle RPA) are successful. This is the case, for example, of the simulation of 
the surface optical spectroscopy, such as reflectivity anisotropy (RA) or differential 
spectroscopy (SDR). 
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On the other hand, new experiments require more complete theories includ- 
ing, for example, spin degree of freedom in order to treat magnetic systems, or 
local field effects in order to describe strong anisotropic systems, or the inclusion 
of semicore and core levels in order to obtain information about core and semicore 
spectroscopies. It is hence important to implement these more complete theories in 
ab initio computational codes in order to describe more realistic condensed matter 
systems. In particular, these implementations are fundamental to be able to treat 
systems with explicit inclusion of surfaces, or isolates sytems. Moreover, new more 
efficient algorithms are essential and the improvement of existing codes is required 
in order to extend the range of numerical simulation applicability to systems with 
a larger size (in terms of number of atoms). 

Theoretical spectroscopy is a successful combination of these quantum theories 
and computer simulation intended to describe the fundamental mechanisms of in- 
teraction between materials and perturbing external fields. 

The present work is an example of what is possible to obtain with the theoretical 
and numerical tools we have just mentioned. In particular, we will discuss problems 
of numerical efficiency and the inclusion of some aspects neglected up to now, such 
as the inclusion of spin variable and semicore levels. 

This manuscript contains different parts: a thread can be drawn from the tech- 
nical development of methods, to simulations of a variety of physical systems. More- 
over, the study of a large variety of complex physical applications helps to point out 
the limits and the advantages of the theories adopted. 

After a brief review of the theoretical background presented in chapter [IJ in chap- 
ter [2] we describe in details the methods used to simulate surface spectroscopies and 
the main experimental techniques. 

From the following chapter, we approach the core of the work developed in this the- 
sis, in paticular in chapter [3] we focus on the dynamical response function %(r, r' ,uj) 
and we present the developement of an Hilbert transform (HT) based method to 
evaluate the independent particle response. The time scaling analysis of the HT- 
method on a model shows that it is convenient for large systems. As an application, 
we studied the crystal local field effects on the optical spectra of the Si(100)-(2x2) 
surface, weakly oxidized. In chapters 0] and El we focused RA and EEL spectra of 
clean and oxidized Si(100) surfaces. Chapter 0] is devoted to the clean surface, for 
which we discuss the spectra calculated on three reconstructions: p(2xl), p(2x2) 
and c(4x2). The oxidation process of this surface is then analyzed in chapter [5l 
where several oxygen adsorption sites are studied through geometric optimization 



Introduction 



3 



and calculation of EEL spectra. 

The following part of this manuscript is devoted to spin polarized systems. The 
limits of DFT and TDDFT-LDA approach are highlighted in the case of BeH, a 
simple molecule with unpaired number of electrons (see chapter [6]). The successful 
calculation of the optical conductivity for bulk iron is then presented in chapter [7] 
and spin resolved electronic properties of this system are provided including the 3s 
and 3p semicore states. 

In the last part of the manuscript, chapter [81 it is summarized the interesting elec- 
tronic and magnetic properties of iron, cobalt and nickel pyrites. These compounds 
are complex spin polarized systems that could have stimulating applications in the 
new field of spin electronics. 
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Part I 

Theory 



Chapter 1 

Theoretical background 



The problem of finding the electronic ground state of a condensed matter system 
is equivalent to solve the fundamental equation of quantum mechanics, i.e. the 
Schrodinger equation for a set of interacting electrons immersed in an external po- 
tential. The Density Functional Theory (DFT) provides a successful tool to treat the 
problem giving a variational reformulation of the equations in terms of the electronic 
density. In this chapter we briefly review this theory. The fundamental theorems 
and its Time Dependent generalization (TDDFT) is reviewed with a presentation 
of the methods used in this thesis for the numerical simulation of realistic systems. 



1.1 The Schrodinger equation for condensed matter sys- 
tems 



The non-relativistic time-indipendent Schrodinger equation of a system constisting 
of N interacting electrons in an external potential generated by M atomic nuclei is 
given by: 

fff(r 1 ,..,r JV )=M(r ll ..,r w ) (1.1) 

where \l/(ri, rjv) represents the wave function of the N-electron many body sys- 
tem. The hamiltonian in eq. (jl.ip is the sum of four operators: 



H 



T e + V ee + V ei + Ti + Vu 

N 



(1.2) 



E 

8=1 



In 1 V- 1 



2 t-r |r, - r jl 



M „ M 

E^kn+E 

a=l a=l 



1 2 l" 



where we assumed atomi units h = m = e = ao = l. The terms in eq. (jl.3p are 
associated to the kinetic energy and the Coulomb interaction of electrons (T e and 
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V ee respectively) and nuclei (Tj and Vu respectively), and to the potential energy 
of the electrons in the field of the nuclei V e i. ^ a (i3) are the atomic numbers of the 
elements involved and are the nuclei distances. 
Several approaches can be adopted in order to solve eq. (|1.3f) . 

First, the pseudopotental approach assumed in this thesis, simplifies the treatment of 
the problem reducing the number of the active electrons just to the valence ones and 
describing for each atom the joint effect of the nucleus and the core electrons with 
a suitable potential. For this reason in the following we will refer to ions instead of 
nuclei in the previous treatment. 

Secondly, ions and electrons masses are extremely differents (M, >> m e ) de- 
termining different time scale motions. Assuming that ions are allowed to move 
adiabatically in the field of the electrons ground state, the problem can be treated 
perturbatively within the Born Oppenheimer approximation. Hence, writing the 
wavefunction as: 

*({r},{i?}) = V{fl}(H)0({i2}) (1.3) 
where {r} = {ri,...,r?v} 
and {R} = {Ri,...,R M } 

it is possible to decouple the hamiltonian separating the ionic and electronic part: 

[T e + V ee + V ei \i> {R} ({r}) = E^{R})^ {Ry ({r}) (1.4) 
[Ti + Vii + E2{{R})]<I>({R}) = Etot<K{R}) (1-5) 

where label i refers to ions and E™ ({R}) in eq. (|1.4p represents the electronic con- 
tribution to the potential energy, i.e. the glue for the nuclei, in fact without this 
attractive term, the system would not be bonded. Conversely in eq. (|1 .5|) ions con- 
tribute to the potential V e i and are seen by the electrons as fixed point charges. 

In conclusion, taking in account the approximations assumed, the Eq. I1.3l reduce 
to the eigenvalues problem of the operator: 

N N 

H = -Y,2 V " + Y. «( r *. r s) - E «('<) (i-6) 

i=l i<j i=l 

where the first term is the kinetic energy of electrons, the third is the external poten- 
tial due to the ions in which the electons are immersed, the second term represents 
the complexity of the problem, it describes the interaction between electrons and 
prevents the decoupling of the equation in N one particle equations. 
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1.2 The variational formulation 

If we consider a time-independent Hamiltonian, as described in the previous section, 
and we assume that periodic boundary conditions are applied, the spectrum of 
eigenvalues and eigenfunctions is discrete. For an arbitrary function ^ with non 
vanishing norm, we can define the quantity: 

E [*] = \\ ' ; 1.7) 

It is then possible to prove that the Schrodinger equation H = E\fy) is equivalent 

to the variational principle 5E = 0. 

In fact, taking the variation of eq. (jl.7j) we obtain: 

8{e m = se [*] (*i*> + e m (<y*i*> + e m (*i<y*) 

= (5¥|H|¥) + (¥|2?|5tf). (1.8) 

and considering that is hermitian we can conclude that: 

6E[V] = O (ff^lff - + - £7|<ftf) = 

O (ff - £ [$] W) = (1.9) 



The importance of the functional defined in eq. (|1.7|) can be seen expanding ^ over 
the wavefunctions {ip n }' 

Xm \ c n\ 2 E n 



=► ff[g]-Eb = E " |C ^ 2(g " (1.10) 

where £Jq = min{E n } is the ground state energy. 
Hence we can conclude that: 

E[V}>E ; $> = a$>o & E[V]= E (1.11) 

i.e. the ground state energy Eq is implicitly defined by the minimization (|1.11|) . 



1.3 Density functional theory 

Density Functional Theory (DFT) is a successful tool largely used in order to study 
ground state properties of condensed matter systems. Within this theory the prob- 
lem of solving the Schrodinger equation for the ground state can be exactly recast 
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into the variational problem of minimizing a functional with rispect to the charge 
density. The complexity of the problem is reduced in principle from having to deal 
with a function of 3N variable to one, the density, that depends only on the 3 spatial 
coordinates. In fact, the key quantity of the theory is the electronic density p(r) 
that, respect to the many-body wavefunction \&(ri, rjy), is a real quantity and 
has an intuitive physical interpretation. A review of the topic can be found in litera- 
ture [IJ El G2 H] , hi this section we will review briefly the most important milestones. 



1.3.1 Hohenberg-Kohn theorem 

The essential role that is played by the charge density in the search for the electronic 
ground state was pointed out for the first time by Hohenberg and Kohn [2]. 
Let us consider a system of N interacting electrons immersed in an external potential 
V ex t with hamiltonian: 

H = H iQt + y ext (1.12) 

in particular H{ nt = T e + V ee and the external potential V^xt is due to the interaction, 
for example, between electrons and ions. Assuming that the ground state is not 
degenerate, the first part of the Hohenberg-Kohn theorem asserts that for every 
density p(r) V-representabkjl], the external potential V ex t is a functional of the charge 
density V ex t = V ex t[p(r)], within an additive constant. 

Let us assume, ad absurdum, that there exists a different potential V^' xt with a 
ground state 1 J/ / corresponding to the same ground state density p(r). If E and E' 
are the respective ground state energies we can write: 

E<{y>'\H\&) = {^'\H'ys>') + {^>'\h - H'ys>') 

= E , + Jp(r)\y ett (r)-V^(r)]dr (1.13) 

A similar equation can be written in case of (ty\H'\*$>), in fact: 

E'<{y\H'\V} = (*|#|tf> + (y\H' - H\V) 

= E - J p(r)[V ext (r) - VUr)]dr (1.14) 

and now adding eq. (|1.13p to eq. f|l . 14j) : 

E + E'<E + E' (1.15) 

A density is V-representable if it is positive defined, normalized to a number N, and such that 
there exists an external potential V(r) for which there is a non-degenerate ground state correspond- 
ing to that density. 
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that is the absurdum. 

Therefore, the theorem establishes the legitimacy of the charge density as the fun- 
damental variable in the electronic problem, demonstrating a one-to-one correspon- 
dence between the density and the external potential p = p[V ex t]- Hence, this 
relation is invertible so the external potential can be viewed as a functional of the 
density V ext = V ext [p\. 

Moreover, since the ground state energy is a function of the external potential 
Eq = EQ[V ex t], it is now possible to write it as a function of the charge density 
(HK functional): 

E HK [p(r)\ = T[p(r)] + E H [p(r)\ + J V cxt p(r)dr (1.16) 
where Eh[p(y)] is the Hartree energy given by: 

E„W)\ = i/*/*<^ (LIT) 

Once that the existance of the HK functional is established the second part of the 
theorem affirms that the minimum of the functional E HK [p(r)] is obtained when the 
charge density p is exactly the ground state density (energy variational principle). 
In conclusion, the total energy of an N interacting particles system can be written 
as a functional of the density. This functional exists, is universal and non depending 
on the form of the external potential, however its analytical form is unknown. 

1.3.2 The Kohn-Sham equations 

The Hohenberg-Kohn theorem provides the theoretical justification to reformulate 
the search for the many-body ground state as a varational problem on the charge 
density. Although the analytical form of the HK functional is unknown, the mini- 
mization procedure lead to a set on N associated differential equations: 



5E HK [p] = 5 



T[p] + V H [p] + J V extP (r)dr - A {^j p(v)dv - N 







m + ^M +v ^ l=x (1 , 8) 

where Vjj is the Hatree potential and A are Lagrange multipliers required by the 
normalization constraint. 

The Kohn and Sham approach is based on the introduction of an auxiliary non 
interacting system of N electrons, having the same density of the real interacting 
system in a suitable external potential. 
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We can write the ground state density of the interacting system expanded on a basis 
of N independent orthonormals orbitals: 



p(r) = ^/^*(r)<Mr) 



11.19) 



where fi represents the occupation factor of the orbital i. 

Hence, the HK functional can be written in terms of the kinetic energy of the non 
interacting system: 



E V*M =T + E H [p] + f V extP (v)dv + E xc [p] 
where To is the kinetic energy of the non-interacting system: 
T [p] = ToM = ^y*0*(r) ~</h(t) 

Eh is the Hartree energy and E xc is the only unknown term: 

E xc [p] = E — T — £"h — V ext - 



dr 



[1.20) 



(1.21) 



(1.22) 



In particular, the E xc term contains the contributions given by the difference in 
the kinetic energy of the interacting and non interacting system, i.e. AT[p] = 
T[p] — Tq[p], the exchange effects (Fermi correlation) and the correlation effects 
(Coulomb correlation). 

Now, we can calculate the minimum of the HK functional (jl.20p for the KS non 
interacting system in a fixed external potential V ex t and with the NxN constraints 
due to the orthonormality of the orbitals: 



N 

m,n 



J mrrri u m,n 







where: 



5p 5 



5(f>* 5ft 6 p Yl 5p 
From eq. (jl.23p we can obtain the following set of equations: 

1. 



V 2 4- V eS 
-V + V DFT 



(pi = K4>i 



'1.23) 



(1.24) 



11.25) 



where H KS = — ^V 2 + Vg?T is the Kohn-Sham Hamiltonian and V§p T is the sum 
of three contributions: 

V$ T = Vh + ^ext + V xc (1.26) 
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the Hartree potential Vh, the external potential V ey ± and the exchange-correlation 
potential given by: 

Vxc = ^ (1.27) 
Now, if we rewrite the density eq. (j!.19[) expanded over the N occupied orbitals: 

N 

p(r)=Y,\Ur)\ 2 (1-28) 

i 

we must solve the set of N one particle equations. Now, if we assume that {4>{\ 
diagonalize the NxN hermitian matrix H KS : 

A m „ = {<p m \H KS \<p n ) (1.29) 

we can write N one-particle equations: 

l~ + Kxt(r) + J y^dr' + V^{v)\ &(r) = A^r) (1.30) 

where A, are now interpreted as the KS energies ef s . 

Since the last two terms of the hamiltonian depend on the eigenvectors throught 
eq. (|1.28p , the eigenvalues and eigenvectors can be determined self consistently. 
The equations (jl.28p and (jl.30p are called Kohn and Sham equations and provide a 
procedure to calculate the total ground-state energy of the system: 

N 

E = E Vc M = £e? 5 - E H [p ] + EM - / Po (v)V xc (r)dr (1.31) 

In conclusion it is worth mentioning that the KS eigenvalues do not have any physical 
meaning, as, for instance, Hartree-Fock eigenvalues that are related to real orbital 
energies via the Koopmans theorem. However, there exists a number of approx- 
imations of the exchange-correlation potential (see sections 11.4.11 11.4.21 and I1.4.3P 
providing good agreements with experimental results in many applications. This 
justifies the practical usefullness of the KS scheme. 



1.4 Methods 



1.4.1 Local density approximation 

Once the Kohn-Sham scheme is defined there still exists the problem of the missing 
analytical representation of the exchange-correlation energy. 
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A commonly used approximation is offered by the Local Density Approximation 
(LDA) [5], which makes DFT practically applicable to a wide variety of systems and 
provides a correct description of systems in which the density varies slowly in space. 
The form of E xc is given by: 

EL DA [p(r)] = J efj>(p(v))p(v)dv (1.32) 

where the local dependence of E xc on the density e X c 9 (p(r)) is given in terms of the 
exchange-correlation energy of the homogeneous electron gas of constant density 
p = p(r). Hence the systems is locally approximated to an homogeneous electrons 
system. The function e X c 9 can be separated in an exchange part: 

and a correlation term € c (p), a function that can be obtained by Quantum Monte 
Carlo simulations (QMC), the most popular form has been given for different den- 
sities by Ceperley and Adler [BJ. 

Moreover the exchange-correlation potential can be written as: 

5E^ DA 



V x L c DA (v) - " 



5n(r) 

= e xc (p(r))+ P (r)^ (1.34) 

The domain of applicability of the LDA has proved to be valid for a large amount of 
systems, even not homogeneous ones. However, its results are not appropriate for the 
case of few electron systems (see chapter [6j. For localized systems self-interaction 
corrections (SIC [7]) are usually used. 

1.4.2 Local spin density approximation 

The Local Spin Density Approximation (LSDA) provides a generalization of the 
LDA to the case of spin polarized calculations. Let us define the spin polarization 
parameter £: 

Q = ft-Pi o<C<l (1-35) 
Pt + Pl 

In the limiting case where £ = 0, p^ = p^ and we will recover the LDA for unpolarized 
systems (U), conversely, if £ = 1 the system is completely spin polarized (P) and it 
is possible to write the following parametrizations (see Ref. [7]): 

e p x (p) = 2 l 'h x , (p) (1.36) 

= ~^(2 4/9 p) (1-37) 
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for the exchange and correlation part, respectively. In the intermediate cases the 
parametrization for e xc is given by: 

e xc (p) = +/(0 - e u xc (p)] (1.38) 

where the smooth interpolation function /(C) is defined by: 

/to - (i+c>4/3 2 : /3 (l _-° a/3 - 2 <^> 

1.4.3 Generalized gradient approximation 

A natural way to improve the LDA in order to account for the inhomogeneities of 
the density is to make a gradient expansion of the exchange-correlation energy with 
respect to the density. In this way e xc results to be dependent on the local derivative 
of the density: 

Eg GA [p h Pi] = J f((>x, Pi, Vp t , V Pl )dv (1.40) 

This is the so called Generalized Gradient Approximation (GGA), often used in 
terms of the Perdew-Burke-Ernzerhof (PBE) [5J E] parametrization. 
The GGA improves the LDA with respect to some applications (molecules or systems 
with strongly inhomogeneous density distribution) but it does not offer a systematic 
advance in the DFT calculation tools. 



1.4.4 Brief review of pseudopotential method 

Pseudopotential approach treats an all-electron variational calculation of ground 
state properties in terms of the only valence wavefunctions immersed in a modified 
potential. In this way, core states, being the most localized and expensives to be 
represented, are not directly included in the calculations: their effect on valence 
electrons is described by a suitable pseudopotential. A review of the topic can be 
found in the literature, ranging from the most influential works [101 [TT1 [12l \13\ [T4"l 
US QUI Ql] to other important but less fundamental ones [T8l fT9| l20l [21] . 
In the following we briefly summarize the most important steps of the method. All- 
electron valence orbitals can be represented as a linear combination of core orbitals 
\tp c ) and a smooth function 

\A) = \^)+J2 a ^ ( L41 ) 

c 

where a cv = (ip c '\(f>^} are coefficients that guarantee the core- valence orthogonality. 
By inverting eq. (jl.4ip with respect to \4>y), it is possible to write valence pseudo 
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wavefunction in terms of all-electron core and valence states. Then, applying the 
Hamiltonian to it is possible to show that they are eigenstates of a modified 
hamiltonian with the same eigenstates of the all-electron wavefunctions: 



The projector defined in eq. (|1.42j) by V = X^ c ( £ f ~~ e c)\'4'c}{'4'c\ is not local. More- 
over, because ((p v \V\<p v ) is positive defined it represents a repulsive and short range 
potential, as it should be to correctly describe core orbitals. 

In the general scheme, norm conserving pseudopotentials are derived from an 
atomic reference state requiring that pseudo and all-electron valence eigenstates 
have the same energies and density outside a chosen core cutoff radius. Normaliza- 
tion of the pseudo orbitals guarantees that they include the same amount of charge 
in the core region. Futhermore pseudo and all-electron logarithmic derivatives agree, 
at the reference energies, beyond the cutoff radius. Finally, norm conservation en- 
sures that the pseudo and all-electron logarithmic derivatives agree also around each 
reference level to first order in the energy. 

In this way a pseudopotential exhibits the same scattering properties as an all- 
electron potential in a neighboorhood of the atomic eigenvalues [22]. This property 
provides a measure of the transferability of the pseudopotential. 

1.5 Time dependent density functional theory 

Density functional theory is a successful tool for a large range of applications, how- 
ever some limits can be underlined. 

First, DFT is a ground state theory and it is not obvious how to generalize the KS 
eigenvalues in order to represent the quasi particle energies^]. Secondly, DFT is a 
theory dealing with stationary states, hence it is not possible to apply it to the case 
of time-dependent hamiltonians. 

Part of those limits are overcomed by the Time Dependent Density Functional 
Theory (TDDFT) that is an exact reformulation of time dependent quantum me- 
chanics where the fundamental variable is the time-dependent electronic density 
p(r, t) instead of the many-body function of the system. The first milestone is the 
Runge-Gross theorem [23] that provide a generalization of the Hohenberg-Kohn 

2 For instance, the direct interpretation of the KS eigenvalues as the quasiparticle energies of the 
system leads to the understimation of the bandgap of semiconductors. 




(1.42) 
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theorem to time dependent densities. The theorem states that there exists a one- 
to-one correspondence between the time dependent external potential W(t) and the 
time dependent density of an evolving system at a fixed initial state: 

W(t)++p(r t t) (1.43) 

If we consider the Hamiltonian describing an N-electrons system given by: 

H(t)=T + V + W(t) (1.44) 

where, beyond the kinetic and the coulombian term (T and V respectively), a time- 
dependent external potential W(t) = ^2iV ex t(r,t) appears that can be expandend 
around an initial time to such as 14xt(r, to) = V^xt(r). 

The time evolution of the system is described by the Schrodinger equation: 

H{t)m = i^(t) (1.45) 

Two time dependent densities p(r,t) and p'(r,t), having a commun initial state 
ip(to) = tpo and influenced by two different external potentials V ext and V^' xt , ex- 
pandable around to and such as V^' xt ^ V ext +c(t), are always different. Hence p(r,t) 
determines the external potential but for a time dependent function c(t). Conversely 
the potential fixes the density but for a time dependent phase: 

Hence, for every time-dependent observable 0{t) that is not depending on time 
derivative or time integral, is a functional of the density: 

m)\o(t)\m) = o[ P ]{t) (i.47) 

From the Runge-Kohn theorem it is straightforward to build the Kohn-Sham scheme 
for the time dependent case (see Refs. [23] and [25]). We can write the action: 

A[j>] = ^ dtm)^ - H(t)m)) (1.48) 

where ip is the many-body wavefunction with initial condition ip(to) = ipo- The 
time-dependent Schrodinger equation corresponds to a stationary point of A[ip\ 
similarly to classical mechanics, where the trajectory is a stationary point of the 
action A = J*^ 1 L(t)dt with L the Lagrangian of the system. The action in eq. (jl.5ip 
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is a functional of the density and has a stationary point corresponding to the correct 
p(r,t), i.e. solving the Euler equations: 

SA[p] 



5p(r,t) 

it is possible to recover the density. 
Similarly to the static case, we can write: 



(1.49) 



A[p] = B[p] - I' dt I fIrp(r.t)V,,,(r.l) (LoO) 
Jt 



where B is a universal functional given by: 



ti g 



B[p] = J t dt{m\i-Q t -T-Vm)) (1-51) 

Now, an auxiliary non interacting system can be associated to the interacting one 
in a similar way than to the Kohn-Sham scheme. The stationary condition can be 
applied to eq. (|1.50p with the condition p(r, t) = ^ \4>i(r, t)\ 2 in order to obtain the 
time dependent KS equations: 



</H(T,t) =i^i{r,t) (1.52) 



~V 2 + V ext (r,t)+ J V(r,r')p(r',t)dr' + V xc (r,t) 

In eq. (|1.52p it is possible to recognize three contributions to the effective potential: 

y cff (r, t) = Vn(r, t) + F ex t(r, t) + V xc (r, t) (1.53) 

the Hartree and the external potential (Vh and V ex t respectively) and the exchange- 
correlation potential defined by the functional derivative: 

KccM) = /^7T (1-54) 

where A xc is the exchange-correlation part of the action (jl.5ip . 



1.6 Linear response 

Within the TDDFT framework we can calculate the linear response of an N particle 
system to an external time dependent perturbation. The response will be related 
the excited states of the system and can be defined as the variation of the density 
with respect to the variation of the time dependent external potential causing the 
perturbation: 

* (My/) = «t^ k - (I ' 55> 
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Similarly, the linear response in the case of the auxiliary non interacting KS system 
can be expressed by: 



X u (r,W 



(1.56) 



5y cfr (r',to IKcff=0 

where the functional derivatives are calculated at the first order in V e xt m eq. (|l-55j> 

and in V e g in eq. (|1.56|) . 

Now, using the following relation: 



Sp Sp 5V e fi o SV eS 
~ X 



we can write: 

SV e s(r,t) 
SV ext (r',t') 

where: 



S(r-r')5(t-t') + 



SV cxt 5V e g SV ext 
S(t - 1") 



SV ( 



cxt 



r — r" 



+ f xc {r,t,r",t'' 



(1.57) 



X (r",t",r',t')dr"dt" 
(1.58) 



Jxc{r,t,r ,t) = — — — — |y oxt =o 



, ) ,/ , 1 1 - 1 (1.59) 

is the exchange-correlation kernel, the quantity that contains the core of the com- 
plexity of the problem. Combining eq. (|1.57|) and (|1.58|) it is possible to write a 
Dyson equation for x and x° : 



X(r,r',w) = x°(r,r',£j) + J dridr 2 x°(r,r 1 ,u) 



1 



ri - r 2 



+ fxc(r i,r 2 , u) 



X(r2,r',a;) 
(1.60) 

The analytical form of the exchange-correlation kernel is unknown, for this reason, 
the solution of this integral equation is not trivial. 

In the case of f xc = the approximation is called the independent particle random 
phase approximation (IP-RPA) that is equivalent to the Hartree theory but with 
the addition of time dependency. In this scheme the density fluctuation at the first 
order is written as: 

p[u) = J x°(r,r',u)V cS (r',u)dr' = J X (r, r', u)V cxt (r', u)dr' (1.61) 

where x° is built using the KS eigenvalues and eigenvectors calculated with an 
approximation for the exchange-correlation potential V xc in the KS hamiltonian. 

The problem of the efficient evaluation of the response function will be discussed 
extensively in chapter [3j where a new method for the calculation of x^ wm be also 
presented. For this reason, we postpone to that chapter the details on the analytical 
form of this quantity. 
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1.6.1 Adiabatic (spin) local density approximation 

The adiabatic local density approximation (ALDA) furnishes a way to compute the 
excitation energies of a system within the TDDFT. Within this approximation the 
exchange-correlation potential defined in eq. (|1.54p is written as: 

V xc (r,t)^ 6 -^M (1.62) 

where the functional derivative is taken respect to the instantaneous density in such 
a way that the exchange-correlation energy depends just on the density at a fixed 
timl By consequence, the exchange-correlation kernel becomes: 

f (rtr' t') - 5V ™( T >Q ~ S (t - t ') SVxcM (I 63) 

f xc {r, t,r,t)- ^ tl) _ d(t t) ^ t/) (1.63) 

and using local density approximation, see eq. (|1.34[) . we can also write: 

f ALDA (r, t, r', = S(t - t')S(v - r') {2^^- + P^§^j (1-64) 

Moreover, if we want to include the spin variable (ALSDA) we obtain the following 
expression: 

f ALDA (r,t,r',t') = 5(t-t')5(r-r>) (1.65) 

( de L J DA (p,Q | de L J DA (p,Q | *e™ DA (p,Q \ 
V dp^ dp± dp^dp x ) 

where the derivation is defined as: 

de(p,Q = de(p,Q | de(p,Q 
dp^ dp dC 

de(p,Q = de(p,Q de(p,() 
dpi dp d( 



1.7 Dielectric function 

The key quantity connecting the theories presented in the previous sections and 
the experimental spectra is the dynamical dielectric function e(r,r',oj). When an 
external perturbing field is applied to the sample, the charge density rearrages and 
an additional potential V m( ± is induced by the polarization of the system. The total 

3 For this reason, in adiabatic LDA memory effects are neglected. 
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potential, or screened potential, is due to the contribution of the external and the 
induced potential: 

Vtot = Vext + V ind (1.66) 
and it can be also written in terms of the dielectric function: 

V to t = yV^rOWrVr' (1.67) 

The dynamical dielectric function can be recovered as: 

e(r, r', oj) = 5(r - r') - J dr"v(r - r')x(r", r', w) (1.68) 

where v(r — r') is the bare coulomb interaction. The dynamical dielectric func- 
tion e{oj) takes in account the rearrangement of the charge density presenting hole 
and charge accumulation due to the perturbation. However, the screened potential 
usually is calculated from the external potential inducing the polarizabilty of the 
system (see eq. (|1.67p ). Hence, the important microscopic quantity is the inverse of 
the dielectric function that can be written as: 

e~ l (r, r', u) = 6(r - r') + J dr"v(r - r") X (r", r', u) (1.69) 

In conclusion, the response function x an d the dynamical dielectric function 
e(r, t',uj) represent the key ingredients for theoretical spectroscopy. In the next 
chapter we will present the connections between e and the three class of experi- 
mental spectroscopies considered in the present thesis: energy loss, reflectivity and 
absorption spectra. 
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Chapter 2 

Surface spectroscopies 



Every real solid is surrounded by surfaces. Moreover the miniaturization of the 
technological devices requires a better understanding of mechanisms at the atom- 
istic scale at which surface effects become important^. 

From the experimental point of view, surface atoms are only visible in sensitive 
techniques or by studying processes involving atoms at the surface (crystal growth, 
adsorption, oxidation, etching, ...). 

Under normal conditions (atmospheric pressure and room temperature) a real sur- 
face of a solid is different from an ideal truncated bulk because of a reordering of 
the surface atomic bonds and because prepared surfaces are normally very reactive 
to atoms and molecules in the environment. From chemisorption to physisorption, 
all kinds of particle adsorption gives rise to an adlayer on the topmost atomic layers 
of the solid. 

Because of this complexity, first principles calculations can be very helpful to 
better understand the physics of such a system. In this section we will briefly 
review the significant experiments and the theoretical tools devoted to describe 
surface physics. 



2.1 Experimental issues 

Spectroscopy is a useful tool to get information about the physical nature or geo- 
metrical reconstruction of surfaces. Many high level experimental technologies has 

In effect if we look at the number of surface atoms (Na(S)) with respect to the bulk (Na(B)) in 
a 1 cm 3 volume cube we can say that surface effects are negligible because: = = 10~ 8 . 

On the contrary in the case of a 100 A length cube we write: ^(v) = = 10 -2 , hence the 
surface signals are not negligible anymore. 



23 



24 



CHAPTER 2. SURFACE SPECTROSCOPIES 



been developed in the last decades in order to create and analyse the surfaces of ma- 
terials, an exhaustive review of the topic can be found in the literature |26[ \27\ 
Here we summarize the highlights to introduce our results presented in the followings 
chapters. 

2.1.1 Preparation and structural properties 

In order to get spectroscopic information, a well defined surface has to be prepared 
on a particular solid, using a special preparation process and under well defined 
external conditions. 

There are several ways to prepare a surface from a crystalline material and they 
can be grouped into three categories: (i) cleavage (limitated to cleavage planes), (ii) 
treatment of imperfect and contaminated surfaces by ion bombardment and thermal 
annealing and (iii) epitaxial growth of a crystal layer by means of evaporation or 
molecular beam epitaxy (MBA). In all cases Ultra High Vacuum (UHV), i.e pressure 
conditions lower than 10~ 8 Pa (10~ 10 torr) are required. 

However, despite the great care in preparing surfaces, irregular deviations from per- 
fect smoothness and purity are always present (steps, terraces or in general surface 
roughness) making real surfaces far from the ideal ones. 

Surface atoms rearrange with respect to the bulk crystal positions because forces 
acting on the on top atoms differ from interactions between atoms inside the vol- 
ume, and as a results this difference can be enhanced depending on the bounding 
behaviour of the material. 

However, the deviation of atom positions from that of an infinite crystal decreases 
with increasing distance from the surface. 

Hence in our theoretical models we will assume with confidence that positions of 
atoms deep inside the bulk are the same as those in an infinite crystal. On the 
contrary, the distortions of the atomic configuration due to the termination of the 
crystal, are important close to the surface. 

In the case of silicon, the main element considered in this work, when a surface is 
created tetrahedral bonds are broken, and a non negligible atomic rearrangement is 
expected to destroy the translational symmetry of an ideal bulk truncated surface. 
Moreover dangling bonds are usually unstable because rebonding lowers the total 
energy pushing surface atoms closer to form pairs (dimers). For this reason we can 
expect that a silicon surface is a good example of a reconstruction process. 
On the contrary, in the case of materials where chemical bonds are less directional 
(as the case of metals), surfaces are created by relaxation of the topmost layers 
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along the direction perpendicular to the surface plane. In this case the changes may 
conserve the translational-symmetry of the bulk. 

The main experimental techniques used in the study of surface structure exploit 
the diffraction of neutral atomic beams or electrons. With Low-Energy Electron 
Diffraction (LEED) the surface periodicity and a reconstruction are observed directly 
via the diffraction pattern, which give an image of the reciprocal lattice. The size 
and shape of the spots contain information about the extention of domains and 
the presence of surface defects. The atomic positions inside the unit cell and the 
relaxation can be studied through the intensity profiles of the diffracted beam, i.e., 
by plotting the measured intensity of each diffraction spot as a function of the 
energy of the incoming electron. The information of the atomic position is obtained 
by comparing the experimental LEED profiles with those obtained by a theoretical 
simulation of the electron diffraction in the crystal, where the atomic positions are 
the input data. An example of LEED patterns is reported in Fig. 12.11 
Reflection High Energy Electron Diffraction (RHEED) is also used, principally to 
monitor the thin film growth. In this technique incident energies of 10 — lOOkeV 
and incident angles of about 3 — 5 degrees are used. 

The structural analysis of LEED is often performed togheter with Auger Electron 
Spectroscopy (AES) to control the chemical composition of the sample. In AES a 
beam of electrons with energies beyond IkeV strikes the surface and the number of 
electron backscattered N(E) is analysed as a function of the energy. The ^ gives 
the signature of the elements present in the sample. 

Moreover we mention the light-ion Rutherford backscattering (RBS), where 
beams of H + or He 2+ ions are used with energies of hundreds of keV (LEIS up 
to 20 keV, MEIS from 20 eV to 200 keV and HEIS to 2 MeV). The ions are scat- 
tered by the nuclei of the crystal following the dynamics of classical Rutherford 
scattering and lose energy along straight trajectories through interaction with the 
electrons. Information on the composition and atomic displacements in the surface 
layers can be obtained thanks to this method detecting the number of ions as a 
function of the energy and the outcoming directions of ions diffused backward. 

Another important technique employed to investigate the structural properties 
of semiconductors is Scanning Tunnel Microscopy (STM). Due to Binnig and Rohrer 
(1982), this technique gives the local density of occupied and empty states integrated 
over a given energy range around the Fermi level. This technique does not need 
UHV. An example of STM image is reported in Fig. 12.11 

From all these techniques we can get structural information about the system 
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Figure 2.1: Examples of experimental images obtained by LEED (left) for a vicinal 
surface, explaining the double spots in the image, and STM (right). In particular in 
a) LEED patterns of a Si(001)-(1 x 2) are reproduced from [31] in b) an STM image 
of the Si(100)-(2xl) surface at 61K, filled (a) and empty (b) states, from [32 1. 

but a comprehensive interpretation of the data must usually be supported by a 
theoretical description which involves the knowledge of the electronic structure [29} 

ED]. 

2.1.2 Electronic properties 

In translationally invariant systems the wave vector k defines a set of good quantum 
numbers for each type of elementary excitation. In the case of an ordered surface 
of a crystal, such a wavevector k, is restricted to two dimensions (parallel to the 
surface) because in the third direction the system is not translationally invariant 
anymore. The Surface Brillouin Zone (SBZ) becomes 2 dimensional and is defined 
as the smallest polygon in the 2D reciprocal space situated symmetrically with 
respect to a given lattice point (the origin) and bounded by points k, satisfying the 
equation: 

k-g = i|g| 2 (2.1) 

where g is a surface reciprocal lattice vector. Figure 12.21 represents three of the 
five SBZ, referring to the ones we considered in the next chapters of this thesis: 
p-rectangular, c-rectangular and square. Further details on surface theory can be 
found in [28]. Among the experimental techniques which allow to inspect directly the 
band structure and the electronic structure in the 2 dimensional BZ we mention the 
most fundamental: Photoemission (PES), Angle-Resolved Photoemission (ARPES) 
and Inverse Photoemission (IPES). 

PES is the most important technique able to give a picture of the Density of 
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Figure 2.2: Surface Brillouin Zone for the three surface reconstructions that we 
considered in the present thesis: a c-rect angular, a p-square and a p-rect angular. 



States (DOS) at the upper atomic planes. The physics behind the PES technique is 
an application of Einstein's photoelectric effect. The sample is exposed to a beam 
of light inducing photoelectric ionization; synchrotron radiation is the ideal isochro- 
matic radiation source. The energies of the emitted photoelectrons are characteristic 
of their original electronic states. For solids, photoelectrons can escape only from 
a depth of the order of nanometers, so that it is the surface layer which is mostly 
analyzed. PES can be performed with X-ray (XES), E~20-150 eV, where it is pos- 
sible to see transitions between surface bands below the edge of the bulk bandgap 
(the small cross section can be improved using grazing angles) |27| . 

ARPES gives information about the k-dispersion of bands and allows to separate 
the contributions from bulk and surface states. The former connect states with the 
same 3D k- vector, the latter involve photoemission processes conserving only the 
component parallel to the surface k||. The detected 2D vector connecting surface 
states and the continuum can be written as: 

k° ut = kf + g Vk ± (2.2) 

where g is a vector of the surface reciprocal lattice. Hence, plotting the electron 
energy as a function of the emission angle 9 by: 



fc||(0) = fc||sinfl = V - " sing (2.3) 

we can say that the peaks in the energy distribution curve represent the initial state 
of the solid labelled by h\. 

Inverse photoemission (IPES) allows to detect the energy of the photon emitted 
when an electron of an external beam of given E and k falls into an empty conduction 
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surface band or an image state. Even optical adsorption is a useful method to study 
the occupied and unoccupied states that in a first approximation can be described 
by the Joint Density of States (JDOS) defined by: 

JDOS{u) = j p(E)6(E F - E)p(E - hcj)9(E + hu - E F )dE (2.4) 

Finally STM is used to describe electronic states of the surfaces by introducing a 
potential difference between the tip and the sample. In this way it is possible to have 
a spatial map of the wave function at different energies for both empty and filled 
states (see Fig. 12. ip . 

In conclusion it is worth mentioning that spectroscopies which study the electronic 
structure are also an indirect test of the surface atomic structure. 

2.1.3 Reflectivity anisotropy experimental spectroscopy 

Optical spectroscopy is an important tool to probe surfaces since they allow for in 
situ, non-destructive and real-time measurements. Moreover, material damage or 
contamination associated with charged particle beams are avoided. 
However, since the light penetration and wavelength are much larger than typical 
surface thicknesses (few A), optical spectroscopy is less sensitive to the surface. 

Nevertheless, a trick can be used in order to resolve the surface signal. This 
is the case of Reflectivity Anisotropy Spectroscopy (RAS) and Surface Differential 
Reflectivity (SDR), optical techinques of great importance for detecting transitions 
between surface states. Surface sensitivity is greatly enhanced with the use of appro- 
priate conditions which enhance the contribution of interband transitions involving 
surface states [33] . 

RAS is defined as the difference between the normalized reflectivities measured at 
normal incidence, for two orthogonal polarizations of light belonging to the surface 
plane: 

Ry ~t~ Rx 

(2.5) 



RAS 



(Ro + ARy) - (R + AR X ) 



2Rq + ARy + AR X 

where Ro is the isotropic Fresnel reflectivity. Since the bulk of a cubic material 
is optically isotropic, any reflectivity anisotropy must be related to the reduced 
symmetry of the surface or to another symmetry breaking perturbation, for instance 
an electric field. In the case of << 1 we can write: 

RAS ~ ARy ~ ARx (2.6) 
Ro 
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Figure 2.3: Experimental RA spectrum of the clean Si(100). We report as an 
example, the results taken from Ref. [33] to show how the choice of the surface can 
affect a RA spectra: Nominal (a) and vicinal (b) surfaces are considered, and it is 
evident how the low energy spectral features are enhanced in the case of the nominal 
surface. 

An example of measured RAS is represented in Fig. 12.31 and geometry scattering is 
shown in Fig. 12. 61 

On the contrary an SDR spectrum is defined by the difference in the reflectivity 
measured on a clean surface before and after passivation (e.g. by adsorbing atoms or 
molecules on the surface). Passivation (oxidation, the case we discuss in chapter [5]), 
removes surface states but does not affect bulk contributions. One hence obtains, 
for the optical response specific to the surface: 

SDR = ARdean ~ AR P ass (2.7) 

-^clean 

In conclusion, it is important to mention that all these techinques can be ap- 
preciably sensitive to the experimental definition of the surface in the sense that in 
some cases (e.g. Si(100) and Si(100):O, as treated in this work: see chapters ??) 
the RA signal is modified because of the presence of steps, terrace or different ori- 
ented domains. The influence of steps on the reflectance spectra has been analysed 
by Jaloviar et al. [35]; hence Shioda and der Weide |36j use highly oriented sur- 
faces (with terraces 1000 times larger than vicinal surfaces) in order to obtain more 
accurate RA profiles. Finally, a comparison of RA spectra obtained by nominal 
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and vicinal surfaces is shown as an example in Fig. 12.31 where we reproduce data 
from Ref. |34j to illustrate how the use of nominal surfaces improves the spectral 
resolution in the low energy region of the spectrum. 



2.1.4 Electron energy loss spectroscopy at surfaces 

A natural complement to optical spectroscopy is Electron Energy Loss Spectroscopy 
(EELS) which, despite some complications in the interpretation of the data, turns 
out to be surface sensitive probe, particularly in High Resolution EELS (HREELS), 
which uses low energy incoming beams. 



Bulk 
Ptasmon 
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Figure 2.4: Schematic repre- 
sentation of all kind of excita- 
tions that can be detected by 
EEL spectroscopy. This figure 
is reproduced from [26] 



In an EEL experiment a material is exposed to a beam of electrons with a defined 
narrow range kinetic energy. Some of the electrons undergo inelastic scattering, 
losing part of their energy and having their paths slightly and randomly deflected 
from the specular direction. The amount of energy loss can be measured via an 
electron spectrometer and interpreted in terms of excitations of the sample. Inelastic 
interactions include phonon excitations, inter and intra-band transitions, plasmon 
excitations, and inner shell ionizations (see Fig. 12. 4[) . 

Usually EELS is performed in transmission, but to be surface sensitive it must 
be applied in a reflection geometry (REELS), see Fig. 12.71 and use relatively low 
incident energies (around 50-100 eV, versus 1 keV for bulk). In HREELS the beam 
is highly monochromatic and the energies of the electrons range up to few eV. 
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2.2 Theoretical surface spectroscopies 

A real surface is a complex physical system whose geometry, i.e. atomic positions, 
is generally unknown and usually involve many degrees of freedom. This makes 
calculations very heavy and can create serious obstacles when fully treating excited 
states. For this reason calculations are usually performed at different levels of so- 
phistication, involving various simplifications and approximations, according to the 
accuracy required and the numerical heaviness. 

We assume to be able to calculate the dielectric tensor of a general bulk system (as 
discussed in the previous chapters) and we are going to describe how to use it in 
order to reproduce and predict surface spectroscopic experiements. 

2.2.1 The slab method 

The description of the crystal termination is solved here using the slab method, i.e. 
representing the surface by means of an atomic slab of suitable thickness (usually 20- 
30 A). Using plane- wave basis sets, the three dimensional periodicity of the system 
can be recovered by considering repeated slabs, separated by a sufficiently large 
region of empty space. 

In Figure 1231 we illustrate the example of a slab inside a supercell and indicate the 
three layers involved (giving the name to the three layers model): a bulk region 
(composed by the inner atoms of the slab), a surface layer with thickness d (top 
layers of the slab), and a vacuum volume. The thickness of the surface layer must 
be smaller than the wavelength A of the light. The bulk properties are assumed to 
be described by an isotropic dielectric function Eb(uj) and the surface is described by 
a frequency dependent dielectric tensor, where the complex diagonal elements are 
defined by e xx (u), e w (w) and e zz (u). 

In practical calculations, the vacuum region is chosen large enough to avoid the 
interaction between the two surfaces of the slab and careful convergence tests have 
to be performed. Figure [231 shows the case of a symmetric slab geometry describing 
an oxidised silicon surface. The slab method is general; also non symmetric slabs 
can be used, even if the case is not treated in this thesis. 

2.2.2 Real— Space slicing technique 

Within the description of the three-layers model, an electron impinging on the sur- 
face feels the potential from this surface layer through its dielectric function e s as 
well as the potential of the bulk region. However, microscopic calculations generally 
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Figure 2.5: Schematic representation of an oxidised silicon surface whithin the su- 
percell approach used in our calculations. 4 oxygen atoms (red circles) and 64 silicon 
atoms (yellow circles) form the symmetric slab. The supercell is the union of the 
16 layer slab and the vacuum region. In figure (a), the bulk, surface and vacuum 
regions of the three-layer model are indicated for the upper half of the slab. On Fig. 
(b) blue line represents an example of a cutoff function used in the real-space slicing 
technique. 

output the dielectric function of the supercell e c . In previous works [37], e s was 
extracted from e c by using the expression: 

J(w) = {N b - 2N s )d l £ b {uj) + 2N s die s {co) (2.8) 

Here d\ is the interlayer spacing, N s (iV&) is the number of layers in each surface 
(bulk) region, and I is the integral of the slab RPA dielectric susceptibility e(uj, z, z') 
over z and z', i.e., I(u) = N c die c . However this approach cannot always guarantee 
perfect cancellation of the bulklike layers in the supercell, and may even lead to 
unphysical negative loss features. A more reliable approach is to extract e s directly 
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using a real-space slicing technique, i.e. by projecting out the response of a defined 
surface layer using a cut off function in real space (for a detailed treatment see 
Ref. [38]). 

This technique is usually reported in terms of the polarizibility of the half slab a hs , 
which, in case of symmetric slabs, is obtained dividing by 2 the full polarizability a: 

A 2 2 

Im[4™£»] = E E \K^\ 2 KEcV ~ E vk - M (2.9) 

k v,c 

where ck are the matrix elements of the momentum operatoi§|. 
We introduce now a cutoff function 6{z) aimed at projecting out the optical transi- 
tions related to a certain selected region of the slab. The function 9{z) is a sum of 
two Heaviside step functions: 

e(z) = H(z-z )-H(z-(z + D c )) (2.10) 

where D c is the thickness of the cutoff function (see Fig. I2.5p . 

The cutoff function is introduced into the calculation of the optical properties 
through the use of a modified matrix element P v k tC \i, defined by: 

^k,ck = -in j drV4(r)^)^ck(r), (2.11) 
and hence the polarizability of the slice is described by the relation: 

Jm[4*c%*(u,)] = E E^ckF^ck^k - E vk - fiu). (2.12) 

k v,c 

In the next paragraphs we will show an application of this technique to analyse the 
layer-by-layer contribution of the electron energy loss spectra. 



2.2.3 Theory of RAS 



By exploiting a reflection geometry and the polarizability of the incident wavevector, 
certain spectroscopies are able to resolve the surface contribution to the optical 
properties of a material. In the case of normal incidence light travelling from a first 
medium with e = £\ to a second medium with e = £2 (see Fig. 12. 6|) the Fresnel 
reflectivity is defined by : 

2 



Rn 



N 2 -N x 


2 


V^2 ~ V^l 


N 2 + N 1 




V^2 + V^l 



(2.13) 



"If the pseudopotential is non local, vm is used instead of the momentum operator 
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Figure 2.6: Schematic illustration of the general SDRS geometry in oblique inci- 
dence (left), and RAS assuming normal incidence (right). The light propagation 
direction is represented by red lines, and we indicate the polarization of light in two 
orthogonal directions belonging to the surface plane. The bulk is assumed to be 
isotropic, while the surface is described by a frequency dependent dielectric tensor 
with eigenvectors parallel to x and y. 



where Ni = yfel is the complex refraction index N = n\ + ini. 
Accordingly, reflectance is defined as a complex number r for which R = \r\ 2 . From 
the experimental point of view, measurements are usually performed with respect 
to this quantity by: 

2 



E rif 



Trine 



(2.14) 



On the contrary, reflectivity is a real number and R G [0, 1] with special cases 
occurring when R = 1 (large difference between N\ and N2), R = (N\ = N2) or 
when Ni for the two media are both real or imaginary. Reflectance is a complex 
quantity linking the electric field amplitudes: 



E = £ e i(/3kx ~ wt) = £ n e^ (Ar:E_rf) = E( ) e i ^ ( - niX ~ ct) e~^ n2X 



oe 



(2.15) 



where oj = vk, and hence k 



-N. From Eq. 12.151 we can write: 



(2.16) 
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from which we deduce that energy decreases exponentially. 

If now we consider the energy density u) and the absorption coefficient a: 

dx 

we obtain: 



-auj 



2lo 



a 
Res 
Ime 



-n 2 



and finally 



C 

nin 2 

m c 



tin 



In the case of an interface with vacuum, the reflectivity is simplified: 



R(u) 



N - 1 


2 


N + 1 





(ni - l) 2 + n 2 



(ni + l) 2 + n 2 



(2.17) 

(2.18) 
(2.19) 
(2.20) 

(2.21) 
(2.22) 



In the case of non normal incidence (see Fig. I2.6P two contributions are distin 
guished: 

/ 2 

cos6> — v e — sin ft 

R 



R„ 



cosft + ye — sin 2 



ecost 



sin 



ecos( 



+ Ve 



sin 



(2.23) 
(2.24) 



related to p and s waves respectively. 

In experiments, the RA spectrum is calculated starting from the knowledge of 
the reflectance and hence the reflectivity by means of Eq. 12. 5[ 
Theoretical models link reflectivity to the dielectric tensor. In the particular case 
where << 1, the relative deviation of the reflectivity with respect to the Fresnel 
contribution is given in terms of the surface and bulk dielectric tensor e s , e s by the 
equation: 



Rq 



— cos#Im 

c 



£6-1 



(2.25) 



representing SDRS formula in s-polarization, or, in the case of normal incidence: 



ARj 

Rq 



4w 

— Im 

c 



Airau(uj) 

£&-l 



(2.26) 



In Eq. 12.261 we used the slab polarizability instead of the surface and bulk dielectric 
function. More details of this theory can be found in reference [33] • 
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2.2.4 Theory of electron energy loss at surfaces 



We use a semiclassical dipole scattering theory that accounts for the long-range 
interaction between the incident electrons and the medium under study [26|, 139] . 
Assuming planar scattering, and taking yz as being the scattering plane (z is the 
surface normal, see Fig. I2.T|) the scattering probability is defined by: 



P(k,k') = A(k,k') Im g(qn, 



UJ 



(2.27) 



where k and k' are the incident and scattered wavevectors respectively. The kine- 
matic factor, A(k, k'): 



A(k,k') 



1 k' 



(eaoir) 2 cosf? k \q?, + q"j_\ 2 



2.; 



mostly contains the information concerning the scattering geometry (see Fig. I2.T[) . 
The angle 9 is the direction of the incident beam with respect to the normal to 
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Figure 2.7: Schematic representation of the reflection geometry for Electron Energy 
Loss Spectroscopy experiments. 



the surface plane, and q\\,q± are the parallel and perpendicular components of the 
transferred momentum q = k — k'. 
The loss function is defined by: 

1 



1 +e cS (qn,ui) 



(2.29) 



and represents the part of Eq. 12.271 involving the approximation of the model and 
the separation between bulk and surface contributions to the dielectric function. 
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If the surface were to be modelled as a semi-infinite truncated bulk, e e g would be 
replaced by and we would obtain the familiar expression of Mills [39J . 

In this work we adopt an anisotropic three-layer model of the surface as derived 
by Selloni and Del Sole [40^ I41j . The surface is modelled as in Fig. I2.8t a semi- 
infinite layer of vacuum, a surface layer of thickness d, represented by the surface 
dielectric tensor e s , and a semi- infinite layer of bulk (dielectric function e^). The 
effective dielectric function is defined by: 



£ e ff(q||,w) = e s (q,w) x 
e 8 (q,u) +e b (g,oj) + A(q, i^e" 2 ^— 
e s (q, u) + e 6 (q, u) - A(q, a^)e~ 2ci l I 



c(q,w) 



(2.30) 



where d is the thickness of the surface, £&(q, ui) and e s (q, ui) are the bulk and surface 
dielectric function and A(q,u>) = £f,(q,u;) — e a (q,a?). 

In particular £ s (q, cj) and the auxiliary function £ aua; (q;w) are written as a function 
of the y, z components of the dielectric tensor: £ s (q, ui) = -y/e Si y(q, £j)£ S)Z (q, w) and 

s- fn ,A - / e s ,y(q,aQ 
£au X (q,Wj - V £a , z (q,a;)- 

Although the dielectric functions appearing in Eq. I2.3UI are fully dependent on qu 
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Figure 2.8: Schematic representation of the constituent parts of the three-layer 
model of the surface. 

and u), such quantities are not easy to calculate, since q and u are not independent. 
Hence we make the approximation of replacing £ s (q, oj) with the optical dielectric 
function £ s (w) ~ lim q ^o £«(q> w). This appears to be a reasonable assumption since 
for most of the experiments modelled in this work, q is rather small. 
Surface dielectric functions are calculated according to Ref. [38] using a cutoff func- 
tion (as discussed in the previous sections) in order to select the number of terminal 
layers contributing to the surface response. 
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Chapter 3 



Efficient calculation of the 
electronic polarizability 

In this section we show the application of an efficient numerical scheme to obtain 
the independent -particle dynamic polarizability matrix x^(r, r', uj), a key quantity 
in modern ab initio excited state calculations. The method has been applied to 
the study of the optical response of a realistic oxidized silicon surface, including 
the effects of crystal local fields. The latter are shown to substantially increase 
the surface optical anisotropy in the energy range below the bulk bandgap. Our 
implementation in a large-scale ab initio computational code allows us to make a 
quantitative study of the CPU time scaling with respect to the system size, and 
demonstrates the real potential of the method for the study of excited states in 
large systems. 

3.1 Motivations 

The recent developments of experimental techniques for the non destructive study 
of solid surfaces call for a simultaneous improvement of the theoretical tools: the 
interpretation and prediction of optical and dielectric properties of surfaces require 
more and more quantitative and reliable ab initio calculations, possibly including 
many-body effects. Such an improvement of the theoretical description can be 
achieved, for example, by lifting some of the usual approximations adopted in the 
calculation of the optical response. However, making less approximations increases 
the computational heaviness, and is only possible if efficient numerical algorithms 
can be adopted. A good example is given by the calculation of the independent- 
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particles dynamical polarizability matrix x ( r > r > which is often required as the 
starting point in Time-Dependent Density Functional Theory (TDDFT) [23] and in 
Many-Body Perturbation Theory-based calculations, such as in the GW [32], or 
GW+Bethe-Salpeter schemes (for a review, see e.g. ref. [3j). Evaluating the full re- 
sponse matrix for realistic, many-atoms systems can be computational challenging, 
since it requires a computational effort growing as the fourth power of the number 
of atoms, and the availability of efficient numerical schemes becomes a key issue. 
Recently, schemes allowing to decouple the sum-over-states and the frequency depen- 
dence have been presented. Miyake and Aryasetiawan [43j and Shishkin and Kresse 
|44] have shown that methods based on the Hilbert transform can substantially re- 
duce the computational cost of frequency-dependent response functions, making it 
comparable to that of the static case. In particular the approach presented in [32] 
has been applied to a linear-muffin-tin-orbital (LMTO) calculation of the spectral 
function of bulk copper, while in [44] , a work focused on the GW implementation 
using the Projector Augmented- Wave method (PAW [IS]), a similar approach is 
used to compute the spectral function of bulk silicon and materials with d electrons 
(GaAs and CdS). Another recent work by D. Foerster [32] is focused on the same 
issue and demonstrates how the use of a basis of local orbitals can reduce the scaling 
of a susceptibility calculation for an N-atom system from iV 4 to N ?J operations for 
each frequency, but at the cost of disk space. 

However, the application of such non-traditional methods to large supercells, such 
as those involved in real surface calculations, have not been presented so far. 
It may be stressed that for a given application the computational burden is deter- 
mined not only by general scaling law, but also by prefactors. In particular, prefac- 
tors determine the crossover where one method becomes more convenient than the 
other. This crossover has not yet been discussed for the Hilbert transform methods. 

In the present work, we demonstrate the application of a scheme -similar to that 
introduced in |43] and |44| - based on the efficient use of the Hilbert transforms, 
by performing the calculation of the optical properties of a realistic, reconstructed 
surface: Si(100)(2x2):O, covered with 1 monolayer (ML) of oxygen. 
We provide a quantitative evaluation of the computational gain for this calculation 
of the full dynamical independent-particle polarizability. The latter is constructed 
from Kohn-Sham eigenvalues and eigenvectors and is then used to compute surface 
optical spectra, including for the first time the local field (LF) effects on Reflectance 
Anisotropy (RAS) and Surface Differential Reflectivity (SDR) spectra of this surface. 
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3.1.1 Local filed effects 

The impact of local fields on surface optical spectra has been a controversial issue 
for decades, specially concerning the so-called intrinsic or bulk-originated effects. 
The latter have been measured, for energies above the bulk bandgap, since the sem- 
inal works by Aspnes and Studna [UJ showing that the normal incidence optical 
reflectivity of natural Si(llO) and Ge(llO) surfaces displays an anisotropy of the 
order of 10~ 3 . The effect was called intrinsic since it is not due to the existence of 
surface states, nor to surface reconstruction. Early model calculations by Mochan 
and Barrera [48j performed for a lattice of polarizable entities and exploiting the 
Clausius-Mossotti relation pointed out that intrinsic anisotropies could be due to LF 
effects. A subsequent work by Del Sole, Mochan and Barrera [37] based on Tight 
Binding (TB) method has shown that the RAS spectra calculated for Si(110):H 
within this semiempirical scheme did not reproduce well the experimental data, de- 
spite the inclusion of surface LF effects. However, more recent calculations based on 
realistic bandstructures (within DFT-LDA with GW-corrected band gap) [19] have 
suggested that intrinsic anisotropies at the bulk critical points for the (almost ideally 
terminated) Si (flO):H surface could arise as a consequence of surface perturbation 
of bulk states, without invoking LF effects. Other tight-binding calculations (see, 
e.g., ref. |50j) suggested the existence of intrinsic surface optical anisotropies not 
due to surface local fields. 

A substantial advance in clarifying the role of local fields has been achieved only 
recently by F. Bechstedt and co-workers, who carried out a calculation of the RAS 
spectra of Si(110):H [51J and monohydride Si(100)(2xl) [52] including self-energy, 
crystal local fields, and excitonic effects from a fully ah initio point of view. In both 
the considered surfaces, which have no surface states within the bulk bandgap, the 
LF were found to cause a slight decrease of the optical reflectivity; however, the ef- 
fect was found to cancel to a large extent in the RAS spectra, being almost identical 
for the two polarizations of the incident light. The situation may be different in the 
case of extrinsic optical anisotropies, i.e. those directly related to surface states and 
surface reconstruction, and appearing below the bulk bandgap. In at least one case 
substantial effects due to LF have been reported [53]. However, further calculations 
for a wider class of surfaces are necessary in order to assess this point more precisely. 
The system we consider here belongs to a widely studied family of surfaces, because 
of their importance in the understanding of silicon-silicon dioxide interfaces in semi- 
conductor technology. Despite the many experimental [Ml ESI [561 EH [581 ESI ED] and 
theoretical [6H [62l [631 EH EHl ES] works appeared in recent years, the debate on the 
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oxidation mechanism of Si(100) is still open. However, the most favorable oxygen 
adsorption sites in the first stages of (room-temperature) oxidation process have 
been identified as the dimer-bridge position, and a bridge position on the backbond 
corresponding to the lower atom of the dimer. This remark is supported by STM 
experiments [60] and by a first-principles molecular dynamics calculation [66]. From 
the theoretical point of view, ground and excited state properties of Si(100)(2x2):O 
at 0.5 and 1ML coverage have been recently studied by some of the authors [67] ; 
however, computational limits prevented till now the inclusion of the local-field ef- 
fects in the ab initio calculation of optical properties. 

In the following paragraphs we brefly summarize the theoretical framework and 
the expression of x^ usually employed in plane-wave based calculations. Then we 
show how the Hilbert transform (HT) technique can be applied, as a generaliza- 
tion of the Kramers-Kronig relations, in order to decouple the sum-over-states and 
the frequency dependence in x^ ■ Moreover an estimation of the accuracy and the 
possible computational gain are presented for a model system. 

3.2 Theoretical framework 

The starting point of our work is a DFT-LDA ground state calculation performed 
with the ABINIT code [68] yielding independent-particle eigenvalues and eigenvec- 
tors within the Kohn-Sham scheme [U [2]. Besides to the occupied ones, empty 
(conduction) states up to an energy of several eV above the Fermi level are obtained 
by means of iterative diagonalization techniques. 

However, in order to study the optical and dielectric response, the level of theory 
must be brought beyond the ground state one, using, e.g., many-body perturbation 
theory or TDDFT [23]. The latter is particularly suited for the study of neutral 
excitations, as those involved in optical reflectivity and electron energy-loss. 

3.2.1 Fundamental ingredients 

Within TDDFT, it is possible to obtain the retarded density-density response func- 
tion x(r, r', uj) from its non-interacting Kohn-Sham counterpart x ( r > r '> w ) through 
a Dyson-like equation: 

x = x (o ) + x (o)^ x (3 _ 1} 
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where the kernel K contains two terms: the Coulomb potential, v c , and the exchange- 
correlation kernel, f xc (r,r',uj). An explicit expression for \ is then given by 



X = x {0) 



,(0) - 1 



1 - (v c + f xc )x W ■ (3.2) 



e M (w) = rim — (3.3) 



Eq.s (|3.ip and (j3.2|) are matrix equations, involving two-points functions such as % 
and ■ I n the present case, working within a plane- waves expansion, XGG'Cli^) 
and XQ G ,(q, u;) are matrices in reciprocal space, and i> c (q + G) = | q + G p * s * ne 
Coulomb potential. The exchange-correlation contribution, f xc , is not exactly known. 
It can be included in an approximate form, e.g. using the LDA functional E] in 
the adiabatic approximation (ALDA), or in a more sophisticated approximation 
such as those described in (BUJ EDI ED E21 E3 El] • In order to compare with optical 
experiments, the macroscopic dielectric function £^(w) must be calculated. The 
latter is denned as: 

1 

where the inverse dielectric function e GG ,(q, u) is linked to the response function 
X by: 

b g,&(.^ u ) = i + ^ c (q + G)xG,G'(q^)- (3-4) 

When only v c is included in the kernel K of eq. (|3.ip exchange and correlation effects 
in the response are neglected, while the use of the correct expression (|3.3p still 
consider the LF effects [75]. Already at this level the calculations can become time 
consuming from the computational point of view when the full x ( r i r ' : U1 ) matrix 
has to be obtained. In complex systems with large unit cells the only tractable way to 
proceed is often to neglect local fields, by assuming that £m(w) is well approximated 
by the average of the microscopic dielectric function: 

eZ LF (u) = lim e ,o (q,^) (3.5) 

q— >0 

This corresponds to neglecting the off-diagonal elements of e in reciprocal space 0. 
When moreover exchange and correlation effects are neglected, (independent quasi- 
particle approximation or IP-RPA) the imaginary part of the macroscopic dielectric 
function takes the simple Ehrenreich and Cohen [76] form: 

lme N M LF {uj) = ^ I < ^l v ^i > -ti-u) (3-6) 

ij 



1 In real space, this corresponds to assume a dependence of ejvr(r, r') only on the difference (r— r') 
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where v is the velocity operator and i, j stand for occupied and unoccupied states 
respectively. The substantial simplification obtained in this case explains why most 
of the calculations of the optical properties of real surfaces are done within the 
independent quasiparticle approach, neglecting local-field effects. On the other 
hand, a fast and efficient scheme to compute the full matrix x^ represents a key 
issue in order to be able to go beyond this approximation, e.g. by including the local 
fields, as we do in the present work. Moreover, an efficient method giving access 
to the full x is °f paramount importance when the screened coulomb interaction 
W GG ' (q) is needed, such as in ab-initio GW calculations. In the following, we hence 
concentrate on the expression of x^ itself, i.e.: 



X (0) (r,r» 



2E^ X " /,)^(r)^(r)V*(r')^(r') 
1 1 



ui - (€j - €i) +ir] oo + (ej - €i) + if] 



(3.7) 



where fi are occupation numbers (0 or 1 in the present case), rj is an infinitesimal 
and the factor 2 is due to the spin degeneracy. Switching to reciprocal space and 
focusing on the case of semiconductors, we make valence (v) and conduction (c) 
bands to appear explicitly, and rewrite this equation as: 
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(e c k - e^k) + iv u + ( £ ck - e«k) + ir l 

(3.8) 

where Oo is the volume of the unitary cell and we have also introduced the notation 
/Sy C k(q+G) to indicate the Fourier transform of </>* k+q (r)(/> c k(r). From the numerical 
point of view the evaluation of these sums for each frequency oo can become very 
heavy. Indeed, for a realistic system the evaluation of eq. (|3.8p involves, for each 
frequency, the summation over a large number of terms, which for a system of 50 
atoms typically is of the order of 10 8 . 



3.2.2 The Hilbert-transform approach 

Since we consider the case of the q — > limit to study optical properties, in the 
following the label q will be omitted to simplify the notation. The generalization 
to the case of finite q is straightforward. Introducing a simplified notation for band 
and k-point indexes, we define a single index of transition t to represent the triplet 
{v, c, k}. In this way, oo t indicates an (always positive) energy difference, (e c ,k — ^,k)- 
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We also introduce the two complex quantities: 

Zi, t = PvcUG)Pouu(G') (3.9) 
Z 2 ,t = -Pcvk(G)pvdc(G') (3.10) 



such that: 



xg» = E 



t 



u> — uJt + if] uj + u t + irj 



(3.11) 



When G = G' (diagonal elements) the Za are real, and Z\ = —Zi- Using 



lim = V ( - ) =F wr<^) (3-12) 



^0+ x±ir] 

one can rewrite the ?? — > + limit of equation ([3. lip as the sum of four terms: 

xfe<«> = < 313 > 

Xgg'M = nr]Tz M {(w-w,) (3.14) 

t 



*&<«> - E^ < 3 ' 15 > 

X G 2 G ,(u;) = i7r^Z 2it «5(a; + a; t ) (3.16) 
i 

R and A label resonant and anti resonant contributions, respectively, and the four 
terms are general complex quantities. In x^ 2 (a;) and x^ 2 (u;) each term Zt con- 
tributes to the function x only at uj = ut, and has no effect elsewhere. By dis- 
cretizing the frequency axis, the sums over t appearing in x R2 an d X A2 can hence be 
performed once and for all, at difference with those labeled by Rl and Al for which 
the sums should be calculated for each uj. Thanks to the linearity of the Hilbert 
transform, defined as 

1 f +co f(x) 
H f(t) = -V \ ^-dx, (3.17) 

71" J-oo X - t 

one can however directly obtain x A1 an d X R1 horn x A2 an d X R2 '- 

X Al = H[ X A2 ] (3.18) 
X R1 = H[ X R2 ] (3.19) 

In such a wajj§, it is possible to recover the complete Xgg'( w ) m the s P ec_ 
tral range of interest from the knowledge of a single sum performed over the poles 



2 In the case of real matrix elements, Zi, n € R, one recovers the Kramers-Kronig relations linking 
real and imaginary parts of the response 
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ojf In other words, one can avoid the explicit summation over t = {c, i>,k} to 
be repeated for each frequency. The present procedure for the calculation of the 
frequency-dependent polarizability matrices is similar to the method of Miyake and 
Aryasetiawan [43] . with the difference that those authors represented ^-functions 
using Gaussians, instead of bare rectangular functions as in our case Q 



3.2.3 Numerical efficiency for a toy system 

Our scheme has been first tested on a model syste nfl in order to check both the 
accuracy and the efficiency of the algorithm. Figure (|3.2p shows the results of 
the test, comparing x^°^(a;) (real part) as obtained in the traditional way (i.e. by 
evaluating expression (|3,8p for several frequencies), and by the Hilbert transform 
(HT) algorithm. The results are practically indistinguishable on the scale of the 
plot. The same figure shows the growth of the required CPU time as a function 
of the number of transitions (number of {v,c, k} triplets). The gain appears to be 
proportional to the system size. The possibility to achieve such a large gain, at 
least in principle and for a simple system, was also noticed in the previous works 
describing efficient algorithms for the calculation of |43|, 144] . 
Alternative approaches for efficient TDDFT calculations have also been suggested. 
In particular, another promising scheme based on a superoperator approach and 
allowing to access TDDFT spectra in a numerically efficient way has been recently 
introduced by Walker and coworkers [77] • This approach is however not designed 
for the calculation of the whole matrix X > contrary to the method studied here. 
In order to know the actual CPU requirements for the calculation of x^°\ an d to 



Figure 3.1: Accuracy test 
for the HT-based algorithm, 
shown for the real part of 
Xqq, (tjj) of a model system 
(see text). The two curves 
turn out to be indistinguish- 
able on the scale of the plot 
(maximum error less than 
0.5% ). 

3 Similarly, Shishkin and Kresse [44] used triangular functions 

4 We considered bulk silicon Kohn-Sham energies, increasing the number of transitions to build 
, and randomly redefining the transition matrix elements 
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explore the possibilities to study complex systems, such as the impurity levels and 
band offsets mentioned in ref . [44J , in practice one has to keep into account the time 
used to compute the matrix elements (numerators in eq. I3.1ip . and the time used to 
perform the Hilbert transforms, which was not explicitly evaluated in previous works. 
In the following, we hence applied our approach, similar in its essence to that used 
in |43j and |44j . to a large system investigating the actual numerical performances 
of the algorithm. As it will be shown below, substantial improvements can actually 
be achieved in such realistic calculations. 
Therefore, in the following section, we use our implementation of the HT scheme 
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Figure 3.2: Computational 
load requested to evaluate 
Xg G ,(w) on a model system, 
function of the number 
of transitions (which scales as 
the size of the system), for 
both the traditional and the 



Size of the system (arb. units) HT-based methods. 

in the ab initio DP code [75] to study a real reconstructed surface: the oxidized 
Si(100)-(2x2), for which we present the first calculation of its optical reflectivity 
spectra (RAS and SDR) with the inclusion of local-field effects. 
Finally, we carefully compare the numerical performance of the DP code with and 
without the use of HTs, and we draw our conclusions. 



3.3 Optical properties of oxidized Si(100)-(2x2) 

The HT method has been implemented into the large scale, plane-waves ab initio 
TDDFT code named DP, developed by the French node of ETSF [79]. As men- 
tioned in the previous paragraphs, we used it to calculate the optical properties of 
Si(100)(2x2):O. For this surface we adopt the equilibrium structure for 1ML cover- 
age shown in figure (|3.4p . which is representative of a situation in which dimer and 
backbond sites are both occupied by an oxygen atom (structure c3 in ref. [801 [HTj 
and see chapter [5] for further considerations). The surface is simulated with a slab 
composed by 6 layers, containing 48 Si and 8 oxygen atoms, in a repeated super- 
cell approach. Our structural results agree well with those of previous calculations 
[62"1 [821 [801 [81]. We use standard norm conserving pseudopotentials of the Hamann 
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type [21], and an energy cut off of 30 Ry, yielding 15000 plane waves in our unit 
cell. Eight special (Monkhorst-Pack, [83]) k-points in the irreducible Brillouin zone 
(IBZ) are used for the self-consistent ground state calculation, while a 7x7 grid is 
used in the evaluation of Xqq'( w )- Kohn-Sham eigenvalues and eigenvectors are 
obtained for all occupied states (120) and for empty states up to 15 eV above the 
highest occupied state (top valence). Optical properties are computed through the 
evaluation of the macroscopic dielectric function with and without the inclusion of 
local field effects. 
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Figure 3.3: Imaginary part of the di- 
agonal components of the slab dielec- 
tric tensor, calculated with (LF) and 
without (IP-RPA) local field effects, for 
the three polarizations: (a) parallel to 
the surface, along the dimer axis; (b) 
parallel to the surface, along the dimer 
chains; (c) perpendicular to the sur- 
face. The spectra presented in this fig- 
ure are not fully converged in the k- 
points sampling. 



Figure [3731 shows the imaginary part of the slab dielectric function as a function 
of the energy. Local-field effects are quite important in the low energy region (0- 
2)eV, enhancing em for light polarized along the direction of the dimers chains (x 
direction, see Figure l37i|) . and suppressing it for light polarized along the dimers axis 
(y direction). This goes in the direction of a better description of the microscopic 
inhomogeneities of the system. In the present case, the extrinsic surface optical 
anisotropy, as defined in the introduction, is hence found to be visibly affected by 
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LF. In the case of the third polarization, i.e. the one perpendicular to the surface 
(not experimentally relevant in the case of normally incident light), local-field effects 
are huge, and introduce a blueshift of the absorption edge as large as 5 eV. This 
can be explained by the strong inhomogeneity of the charge distribution in passing 
from the slab to the vacuum, leading to a classical depolarization effect. Similar 
behaviors have been found for example in GaAs/AlAs superlattices |84j . in graphite 
[85] and nanowires [HH [87] . 

Starting from the slab dielectric function, we computed Reflectance Anisotropy 
(RAS) and Surface Differential Reflectivity (SDR) Spectra [88], with and without 
inclusion of LF effects. We used theoretical models (see chapter [2] or [53]) linking 
the RAS and SDR spectra to the dielectric functions evaluated for the bulk crystal 
(s&) and for the slab (en) through the relation: 



where e xx and e yy are the diagonal components of the surface dielectric tensor. We 
show our results for RAS and SDR in figures 14.91 and 13.61 respectively. 



A _ 4Wj ~ £yy(u) - e xx (m) 
R c ' I e b (u) 



(3.20) 




Figure 3.4: Surface struc- 
ture of Si(100)(2x2):O at 1ML 
coverage with oxidation of Si 
dimers and backbonds, oxy- 
gen atoms are depicted in 
dark gray (red), while light 
gray (yellow) circles represent 
bulk and surface Si atoms. 



along the x direction, (a): top 
view of the surface (xy plane), 
with the surface unit cell; (b): 
lateral view of the half slab 
(yz plane). 



Dimer chains are oriented 



We first discuss the case of RAS. The effects of local fields on the imaginary part 
of the dielectric tensor are most evident in the low-energy region of the spectrum 
(below 2 eV), as shown in the inset of Figures |3.3r and 13.3b . In particular, LF are 
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Figure 3.5: Calculated RAS spectrum of Si(100)(2x2):O at convergence, for the 
structural model shown in figure I3T41 Results including (LF) or neglecting (IP-RPA) 
the local field effects are very similar, except for the region between 0.8 and 1.8 eV, 
where the LF effects strongly enhance the RAS signal. The energy scale has been 
shifted by 0.6 eV to compensate for the neglect of self-energy effects. 



found to enhance and sharpen the strong £2 peak at about 1.2eV for light polarized 
along the dimer chains (fig. 13.3b ). and to reduce the first three peaks for light 
polarized along the dimer axis. As a result, LF induce a strong enhancement in the 
surface optical anisotropy (of the order of 100%) in the region between 0.8eV and 
2eV, as displayed in fig. 14.91 This low-energy region (below the direct gap of bulk 
Si) corresponds to surface-localized states, which are expected to carry the surface 
anisotropy. The fact that LF evidence this anisotropy is consistent with the fact 
that dimer chains realize a structure which is geometrically strongly inhomogeneous 
in the direction perpendicular to the dimer chains (see fig. 13. 4p . At higher energies 
(above 2eV) bulk contributions dominate £2, and the resulting RAS is mainly due 
to surface perturbed bulk states. The latter appear to be less affected by local fields 
than the true surface states, and lead to a RAS spectrum which, above 2.0eV, is 
almost insensitive to the inclusion of local field effects. This picture is confirmed 
by the analysis of SDR results. The latter are in fact calculated for unpolarized 
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Figure 3.6: Calculated SDR spectrum (unpolarized light) of Si(100)(2x2):O, for 
the structural model shown in figure I3.41 Results including (LF) or not (IP-RPA) 
the local fields are almost indistinguishable, showing that the effects visible in the 
low-energy part of figure I3.3I are canceling each other in the SDR spectrum. The 
same energy shift as in figure I3~9l has been applied. In polarized SDR the local field 
effects would be of the same size as for the RAS spectra. 



light, i.e. by averaging e xx and e yy . Since LF enlarge e yy and reduce e xx , their 
effects almost completely cancel out when the average is taken. Our calculated 
(unpolarized) SDR spectrum, displayed in fig. 13.61 appears in fact to be very little 
affected by the local fields, in the whole energy range between and 6 eV. However, 
if a polarized SDR spectrum is computed, then local fields are found to influence 
the low-energy region (< 2 eV), in a way which is very similar to the behavior of 
the RAS. 

Unfortunately, it is not possible to perform here a comprehensive comparison 
with RAS and SDR experimental data, since this would require the calculation of 
several possible reconstruction and geometries. In fact, the oxidation mechanism 
of Si(100) has been shown to be exceedingly complex, with different mechanisms 
playing their role depending on the oxidation temperature: a barrierless oxidation 
of the first Si layer [66] , or an "active oxidation" involving etching of the surface and 
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penetration of oxygen in a layer-by-layer manner at higher temperature [891 l90|l9T]. 
Recently, the Si(100)(2xl):O surface optical anisotropy has been shown to be sen- 
sitive to the structural details of the oxygen adsorption by ab-initio calculations of 
the atomic geometries and optical response of a large number of Si(100):O struc- 
tures [92l E3]. An highly structured potential energy surface has been found, with 
minima at the backbonds of the "down" atoms in Si-Si dimers [U3]. Moreover, an 
appreciable amount of disorder is probably present after oxidation of the first Si 
monolayer, and the local strain induced by oxygen adsorption is expected to have a 
sizable impact on the optical anisotropy spectra [92J. 

However, our findings for LF effects in the single case studied here suggest an 
important general remark about surface optical spectroscopies. In fact, it is well 
known that, due to the large penetration depth of visible and UV photons, the 
surface-specific optical reflectivity signal is very small with respect to the bulk con- 
tribution. For materials with an isotropic bulk, the RAS spectroscopy has indeed 
been developed in order to extract the surface signal, by exploiting its anisotropy. A 
correct evaluation of the latter has hence the highest priority in theoretical calcula- 
tions of surface optical spectra. The fact that crystal local fields are potentially able 
to alter significantly the surface optical anisotropy, at least below the bulk bandgap, 
should hence be kept in mind, particularly when the anisotropy of electronic states 
is associated with a large structural anisotropy at the surface, such as in the case of 
dimer chains on Si(100)(2xl). 



3.4 Computational scaling and performances 

In this section, we present a quantitative analysis of the numerical performance of the 
HT-based approach, as implemented in the large scale code DP [75], with respect to 
the traditional approach. Several calculations have been done by varying the three 
main convergence parameters: (i) the number of valence-conduction transitions 
(Nt = Nk x N v x N c ); (ii) the number of frequency intervals considered in the 
spectrum, i.e. the spectral resolution (number of frequencies, N u ); (iii) the number 
of plane-waves considered in the response matrix (Nq). Optical properties usually 
converge at an Nq value which can be substantially smaller than the total number 
of plane waves, N g , used to describe the wavefunctions. 

The calculation of expected to scale, in the case of the reference approach, 
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T ref = N k [aN v N c N g logN g + PN v N B N%iN u ] + A (3.21) 

where N k is the number of k-points, and a and /3 are prefactors which are indepen- 
dent on N v , N c , Ng, Nk, N g , and N u . The first term in equation (|3.2ip is due to the 
evaluation of the numerators Z n in equations (|3.14|) and (|3.16|) by using FFT, and 
is present both in the reference and Hilbert approach. The second term stems from 
to the evaluation of equation (|3.1ip in the traditional way. The remaining term A 
keeps into account residual parts of the calculation, as the matrix inversions, which 
contribute much less to the CPU time than the first two terms. 

The expected scaling in the case of the Hilbert-based scheme is instead: 

T new = N k [aN v N c N g logN g + P'N V N C N G ] + jN G N% + A' (3.22) 

In this case, the second term does not contain the factor N& anymore, and its 
prefactor becomes /?', due to the calculation of \ A2 an d X R2 (equations 13.141 and 
I3.16p . The calculation becomes, in this sense, comparable to a static one. 

However, the actual evaluation of the Hilbert transforms (equations 13.181 and 
I3.19P introduces a new term scaling as N G N%. Due to the small prefactor 7, the 
latter term can often be neglected (see, e.g., figure 3 of reference |43|). In the present 
work, we found that the CPU time spent inside the Hilbert transform itself can be 
made negligible by an optimized algorithm^]. 

Considering an N-atoms unit cell, the number of transitions Nt is clearly the 
parameter growing fastest with the system size, since it is proportional to N 2 . About 
22000 transitions per k-point, corresponding to the inclusion of about 200 empty 
bands, are requested to converge the dielectric tensor of the Si(100)(2x2):O slab 
up to 12 eV. The number of frequency intervals, N u , is instead independent on 
N, but it grows linearly with the required spectral resolution. In the present case, 
300 frequencies have been necessary in order to achieve a 40 meV resolution over 
a spectral range of 12 eV. Finally, Ng, i.e. the size of x {0) in reciprocal space, 
depends on the requested real-space resolution needed in the description of the 
induced density variations. This means that larger Ng will be necessary to describe 
systems with smaller interatomic distances, or with larger polarizability. The real- 
space resolution is independent on the system size; however, for a fixed resolution 

5 A r tA f Q leads to the N%t scaling mentioned in the introductory section. 

6 Exploiting the fact that the principal value numerical integration routine has to be called JVq 
times, always on the same energy intervals, a substantial gain could be achieved by tabulating the 
(about 10 6 ) required values of the complex logarithm once and for all at the beginning of the double 
loop over the number of G vectors in the construction of x' '- 
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Figure 3.7: Quantitative 
study of the computational 
load required to evaluate 
the matrix for the 56- 

atoms slab representing a 
Si(100)(2x2):O surface. The 
effects of the three main 
parameters determining the 
numerical convergence of 
the theoretical spectra are 
studied separately, (a): num- 
ber of valence-conduction 
transitions, determined by 
the energy cutoff on the 
empty (conduction) bands 
included in eq. (|3.8p . the 
number of occupied bands 
being fixed; (b): number of 
frequency intervals taken on 
the u axis, determined by the 
requested spectral resolution; 
(c): size of the matrix 
in reciprocal space, roughly 
proportional to the system 
size (for a fixed real-space 
resolution). 



200 300 
Ng {l matrix size) 



Nq, will grow linearly with the volume of the unit cell in direct space. A^g is hence 
proportional to the number of atoms (in a bulk system) or to the volume of the 
supercell (for a finite or semiinfinite system). In our slab calculation, converging 
the spectra with local field effects requires to consider at least 113x113 matrices 
(incidentally, we stress that the IP-RPA spectra, requiring just the G = 0, G' = 
matrix element of x \ do not depend on iVc )• 
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The HT algorithm turns out to be clearly advantageous with respect to Nt and 
Ng, as shown in the first and last panels of figure 1X71 Concerning N w , despite an 
unfavorable scaling in the limit of infinite spectral resolution (the CPU time grows 
quadratically with the number of frequency intervals), one must notice that, due to 
the small prefactor 7, the HT method is largely convenient in whole range of interest 
(N u ~ 10 3 ). 

3.5 Conclusions 

Two classes of conclusions can be drawn from the presented results. First, the 
physics: the study of Si(100)(2x2):O has shown that local field effects, although 
playing only a minor role on the surface optical properties above the bulk bandgap, 
are able to enhance substantially the surface optical anisotropy in the low-energy 
end of the spectra. A similar effect can be expected for other surfaces, when the 
anisotropy of electronic states is associated with a large structural anisotropy, such 
as in the case of the dimer chains on Si(100)(2x2). Moreover, large local field 
effects are found for light polarized normally to the surface. Second, the numerics: 
the computational gain achievable by using the Hilbert transform-based algorithm 
has been shown to be substantial, both in a model system and in a real, physical 
application. A successful implementation of the Hilbert transform method in the 
large-scale plane-waves ab initio computer code DP [78] allows us to locate the 
crossover (starting from which the Hilbert transform algorithm becomes convenient) 
already at medium size systems (less than 50 atoms). 
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Chapter 4 

The clean Si(lOO) surface 



Surfaces are complex physical systems that are very important from both the fun- 
damental and technological point of view. One of the most important surfaces of all 
is the (100) surface of silicon, as it appears in most electronic devices. Nevertheless, 
in spite of its many applications, a completely unambiguous structural and physical 
description of Si(100) is lacking. In this chapter we summarize our results about first 
principles calculations of electron energy loss spectra of the Si(100) clean surface. 

4.1 Si(100): which reconstruction? 

Due to its enormous technological importance, the Si(100) surface has been the sub- 
ject of a wide range of experimental and theoretical studies spanning several decades. 
In fact, quality publications continue to appear regarding the atomic structure and 
electronic properties of the clean surface. Following early LEED experiments |94j . 
it was understood that Si(100) forms a p(2 x 1) reconstruction. The classic ex- 
planation of the LEED observation is that the surface is composed of rows of Si 
dimers separated by trenches (Fig. 14.2]) . as confirmed by various scanning tunnelling 
microscopy studies [951 IHS] ■ Although some quantum chemistry studies have found 
that a symmetric dimer structure (causing a metallic surface) forms the global min- 
imum |97j . several total energy calculations based on density functional theory |98j 
have found that dimer buckling induces a small energy gain, such that the dimers 
adopt an asymmetric configuration and the surface remains semiconducting [99J. 

Three distinct structures have been proposed for the Si(100) surface: the p(2 x 1), 
whereby all dimers are buckled the same way (see Fig. 14. 21( a)): the p{2 x 2) structure, 
where alternating dimers in a row are buckled in opposite directions, and adjacent 
rows are buckled in phase (see Fig. 14.2( b)): the c(4 x 2) phase, being the same as the 
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p{2 x 2) but with adjacent rows buckled out of phase (see Fig. 14.2( c)). Total energy 
calculations have found [98J that the p{2 x 1) reconstruction is prohibitively higher in 
energy than the other two (at zero K), and that the c(4 x 2) is only slightly favoured 
over the p(2 x 2). Which reconstruction occurs on the surface depends critically on 
the temperature. LEED studies have shown that an order-disorder phase transition 
occurs at about 200K [1U(J|, I1U1] . Below this critical temperature, a c(4 x 2) phase 
is generally observed; above it, a p(2 x 1) periodicity is seen. Direct observation of 
the surface structure with STM is complicated by two factors however. Firstly, it 
is well established now that the experimental measurement itself can influence the 
result, and drive c(4x 2)—>p(2 x 2) phase transitions \102\ F103j. Charge injection, or 
electric fields induced by the STM tip, can cause dimers to flip, according to various 
experimental |104t 1103] and theoretical works [105} 1106] . Nevertheless, the general 
consensus is that the c(4 x 2) reconstruction is the more stable structure below the 
critical temperature [W3\ I1U7] . 

Above 200K, STM images appear to show a symmetric dimer configuration. 
However, at these temperatures the dimer rocking mode is activated, and hence it 
is believed that the observed symmetric p{2 x 1) structure is merely a time average 
of the thermal flip-flop motion of the buckled dimers. Based on molecular dynamics 
simulation of the dimer motion, it was suggested that the surface consists of simul- 
taneous local presence of asymmetric dimers, and of instantaneously flat symmetric 
dimers [108] . More recent studies have suggested that the dimers remain short-range 
correlated (see refs. 32-35 in Ref. [109J). In particular, a two photon photoemission 
(2PPE) study found minor difference between the surface band dispersion at 90K 
and at room temperature [1 10] . 

Theoretical simulations of the reflectance anisotropy (RA) spectra confirm that 
the p{2 x 1) reconstruction does not reproduce correctly the experimental line- 
shape 11121 1113] . On the other hand, both c(4 x 2) and p(2 x 2) structures 
are quite similar to the experiment, with the c(4 x 2) yielding a slightly better 
agreement |111| . in particular predicting the observed SDR structure below 1 eV. 
Simulation of the surface differential reflectance also favours the c(4x2) surface [114J. 

In addition to optical techniques such as RAS or SDR, electron energy loss spec- 
troscopy (EELS) in the reflection geometry (REELS or RELS) offers an enhanced 
surface sensitivity and easy access to a wide energy range (see chapter [2]) . Although 
the majority of literature considering REELS of Si(100) has focused on vibrational 
properties [39], several studies have examined the nature of electronic states at the 
clean Si(100) surface [115] . Indirect information about surface states was derived 
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from related studies looking at the changes in the REEL spectrum following oxida- 
tion [1161 1117] 1118] . High resolution (HREELS) measurements were carried out by 
Farrell et al. [1 19] and Gavioli et al. [120] . In the latter work, tight binding calcula- 
tions were performed on the p(2 x 1) symmetric and asymmetric dimer models, and 
suggested that a mixture of the two structures was necessary to explain the room 
temperature EEL spectra. However, the c(4 x 2) structure was not considered in 
that work, and therefore some of the conclusions reached are not complete. 

In the following sections, we present a computational study of high resolution 
EELS for the different reconstructions of Si(100): p(2 x 1), p(2 x 2) and c(4 x 2) We 
consider the energy range that probes the excitation of interband, i.e., about 0-6 eV, 
and hence connect the experimental observation directly with the atomic structure 
and microscopic electronic response. 

4.2 First principle scheme 

We use density- functional theory within the local density approxiation (DFT-LDA) , 
within a plane-wave and pseudopotential framework. The ABINIT |121|, I122j and 
quantum-espresso/PWSCF |123j codes were used for computing the relaxed atomic 
structures, electronic bandstructures, and Kohn-Sham eigenvalues and eigenvectors 
required for the evaluation of the optical properties. However, in order to be consis- 
tent and since we found only minor differences between spectra computed with the 
two codes, we report only the final results obtained using the output of PWSCF. 
For all our calculations we used standard norm conserving pseudopotentials of the 
Hamann type [2T] for both silicon and oxygen generated with FHI98PP package 
|124j within DFT-LDA (Perdew-Zunger parametrization [7]) framework. 

4.3 Geometric structure 

Even if a 7.5 Ha kinetic energy cutoff is rasonable in order to treat the case of a 
clean silicon surface, we used a 15 Ha cutoff throughout, as was previously found to 
be sufficient for oxidized Si(100) with the same pseudopotential [67]. The effect of 
the two cutoff is shown in Fig. 14.11 where it is clear that the minimum of the curve 
is not appreciably different. 

Therefore we used the theoretical lattice constant of 5.393 A, as determined at 
15 Ha A standard repeated-slab and supercell approach was adopted in order to 
model the surface structure. 
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Figure 4.1: Optimization of the lattice parameter of bulk silicon at 7.5 Ha and 
15 Ha cutoff energy. 



We use relatively thick slabs (16 atomic layers) separated by 8 layers of vacuum 
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Figure 4.2: Ball and stick model of p(2 x 1), p(2 x 2) and c(4 x 2) surface recon- 
structions. Large circles indicate "up" silicon dimer atoms. Unit cells are indicated 
by shaded regions. 

(about 10 A). During the geometry optimization, the central four layers were fixed at 
the bulk positions and structures were relaxed until the cartesian force components 
were less than 20 meV/A. Our obtained structural parameters are similar to those 
obtained previously for this surface (see for example Ref. [Ill] ) such as a dimer 
buckling of 0.755 A and a dimer length of 2.33 A. 
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4.4 Computational aspect 

Optical spectra and energy loss spectra were calculated using the YAMBO code 
[125j . taking Kohn-Sham eigenvalues and eigenvectors from a non-self-consistent 
run of the PWSCF code [123J. We carried out convergence tests on the optical and 
energy loss properties with respect to the number of bands and the number of k- 
points. For the calculation of energy loss we included up to 320 unoccupied bands, 
corresponding to minimum e—h transition energy of ~ 13.5 eV. Further convergence 



SiMDD! p(2x2} 36k Ibach IOOcV 




Figure 4.3: EEL spectrum of the clean Si(100)-p(2x2) surface calculated for in- 
cident energy Ej=100eV and 6 = 42. The height of the surface plasmon peak is 
mostly affected by the convergence respect to the number of bands. 



tests were carried out regarding k-point sampling. For the most converged calcula- 
tions, we used roughly equivalent density sets for the three reconstructions: for the 
c(4 x 2), 72 k-points in the irreducible Brillouin zone (IBZ), which is equivalent to 
1152 points in the (1 x 1) BZ (for both clean and oxidized cells, see next chapter); 
for the p(2 x 2), 64 k-points in the IBZ (equivalent to 1024 points in the (1 x 1) BZ); 
for the p(2 x 1), 100 k-points in the IBZ (equivalent to 800 points in the (1 x 1) BZ). 
All spectra reported in this work were obtained using the non-interacting particle 
(RPA) level of theory. Many-body effects, including local field and excitonic effects, 
were compensated for by applying a scissors operator of +0.5 eV to the unoccupied 
states, following the recipe of Del Sole and Girlanda [37]. In this way we account for 
the well-known underestimation of the DFT band-gap, and partially include self- 
energy and excitonic shifts in energy. The value of +0.5 eV was determined in other 
works on Si(100) as giving best agreement with the experimental RAS [67] . 
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In the following subsection we present some discussion of major technical tools used 
in order to compute REELS. 

4.4.1 Broadening 

Inspection of Eq . 1X251 reveals that the cross section goes as 1/q 3 as u — > 0. Hence any 
features appearing in the loss function Im g(q» , oS) at low energy can be dramatically 
enhanced by the kinematic factor. From the computational point of view this means 
that unphysical features may appear close to the origin if a Lorentzian broadening 
is used when calculating the surface dielectric function of a system with a small 
band gap. We adopt therefore a tiny Lorentzian broadening of 6l = 0.006 eV when 
calculating the dielectric function, and afterwards convolute the loss spectra with a 
Gaussian (FWHM = 0.3 eV) to approach the experimental resolution. 

4.4.2 Slicing methods 

A crucial adjustable parameter present within the three-layer model of energy loss is 
the thickness d of the surface layer, as used in Eq. 12.301 According to the model, an 
electron impinging on the surface feels the potential from this surface layer through 
its dielectric function e s as well as that of the bulk layer e\,. Obtaining e s from the 
dielectric function of the slab or supercell e c is not trivial. However early efforts 
used a simple subtraction of the computed bulk dielectric function |37j : 

I(u) = (N b - 2N s )d l e b (co) + 2N s d l e s (co) (4.1) 

where I is the integral of the slab RPA dielectric susceptibility e(w, z, z') over z and z', 
d[ is the interlayer spacing, N s is the number of layers at each surface with dielectric 
function e s (oj) and N& is the number of inner layers with bulk dielectric constant 
Sb(oj). Unfortunatly this approach cannot always guarantee perfect cancellation of 
the bulklike layers in the supercell, and may lead to unphysical negative loss features. 
A better approach is to extract e s directly using a "cut off' function as described in 
chapter [2] or in Ref. [38] , Nevertheless, the choice of d remains somewhat arbitrary. 
In this work we define the lower bound of the surface layer corresponding to the 
actual penetration depth of the electron for the chosen incident kinetic energy, with 
the upper bound defined by the maximum extent of the slab charge density or a 
typical surface state wavefunction (see Fig. 1-4. 5|) . This turns out to be roughly one 
atomic layer. We then checked that varying d by an atomic layer did not change 
the results too much. The dependence of the REEL spectra on the parameter d is 
illustrated in Fig. 14.51 
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.4: Effect of the Lorenzian broadening (top), effect of the gaussian con- 
bottom) on the EEL spectra of a Si(100)-c(4x2) surface (Ej=7 eV and 



4.4.3 Detector integration 



In experimental EEL spectroscopy the detector has a finite acceptance angle for 
collecting the scattered electrons. Hence, in order to improve the quantitative es- 
timation of the EEL spectra at the surface, we implemented a method to perform 
a numerical integration over the circular detector. In particular we account better 
for the finite size of the detector window through a random sampling of the circular 
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Figure 4.5: Schematic diagram of the Si(100)-c(4 x 2) slab and calculated layer 
averaged charge density (left); atomic layers are marked with horizontal lines. De- 
pendence of REELS spectrum on the chosen surface layer (right). 



area and in the computation of the following integral: 

A(k, k')lmg(qii,uj)dQ 



(4.2) 



We implemented the method in the Yambo code [126J , the code we used for all fur- 
ther calculations. In Fig. I4.7I an illustration of the importance of numerical detector 
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Figure 4.6: Integration geom- 
etry. 



integration scheme with regard to relative intensity of peaks is shown. A frequently 
used technique is to approximate the integral by an averaged value of the transferred 
momentum q. 
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Energy (&V) 

Figure 4.7: Dependence of REEL spectra on detector integration method: numer- 
ical Monte Carlo, mean value scheme, and single point sampling. Example shown 
for c(4 x 2) surface, E = 40 eV; O = 60°; 6> det = 1°.. 




Figure 4.8: Dependence of 
REEL spectra on the detector 
size. 
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A comparaison between the two methods is shown in figure 14.71 where relative in- 
tensity of the low energy peaks is corrected with converged numerical integration. 
We also add a figure showing the effect of the detector size on the EEL spectrum. 

4.5 Calculated spectra for Si(100) 

A review of the current understanding of the structure of the clean Si(100) surface 
was given in the introduction and a schematic diagram of the three basic surface 
reconstructions (the p{2 x 1), p(2 x 2) and c(4 x 2)) is given in Fig. 14.21 Throughout 
this work we will refer to the [Oil] direction as x and the [Oil] direction as y, with 
[100] being the surface normal z. Here we present our results concerning optical and 
electronic spectroscopic study of the clean Si(100) surface. 
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4.5.1 Reflectivity anisotropy spectra 

Several theoretical studies of the RA spectra for the Si(100) surface have previously 
been carried out, including tight binding calculations [127] . discrete-dipole models 
\12S\ I129| . ab initio calculations at the independent particle level |130[ I13H 1132} 
II 141 [671 11121 1113] as well as more recent studies including many-body effects [133] . 
Generally it was found that the best agreement with the experimental RA data [36] 
is obtained when the c(4 x 2) or p{2 x 2) models are used in the calculations, while 
the p(2 x 1) gives poor agreement. 

Since we will use the same optical dielectric functions when computing the energy 
loss spectra, we show for completeness (in Fig. 14. 9p the results of our own supercell 
calculations for the c(4 x 2), p(2 x 2) and p(2 x 1) reconstructions at the independent 
particle level. The reflectance anistropy (see chapter ED is defined as 

RAS = ^_^, (4.3) 

where ARi/R (i = x,y) is the normalized reflectivity (i.e., relative to the Fresnel 
reflectivity) . 

As expected, we find that both c(4 x 2) and p(2 x 2) spectra yield a good agreement 
with the experimental data. 

Thep(2x 1) reproduces the low energy peak at 1.5 eV rather well, but the comparison 
worsens at higher energy. Unfortunately, no experimental data is available for the 
RAS of Si(100) in the near-IR range, and is hence limited to >1.1 eV. 

4.5.2 Calculated REELS spectra at E = 40 eV 

We now contrast the RAS results with the theoretical simulation of the HREELS 
experiment of Farrell et al. [119] with incident energy of Eq = 40 eV, which roughly 
covers the same spectral range as the available RAS data. Experimental results 
refers to specular geometry with incident angle of 9 = 60°. 

There exists a one-to-one correspondence between RAS and HREELS that can 
be justified with theoretical arguemnts when we consider a spectral range below E\. 
In fact it can be written that RAS is proportional to Ime s when Ime^ is close to 
zero. On the other hand: 



\mg{q, to) ~ Im 



1 



l + e s 



lm£ ° (4.4) 



'1 + Ree s ) 2 + line? 



giving a direct proportionality of REEL respect to Ime s . This correspondence is 
evident in Fig. 14.141 where HREEL peaks and Ime s are represented. The experi- 
mental one-to-one correspondance between RAS and HREELS was also previously 
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Figure 4.9: RAS spectra of the c(4 x 2), p(2 x 2) and p(2 x 1) reconstructions 
of clean Si(100), compared with the experimental spectrum of the nominal surface 
(scaled by a factor of 5). Experimental data are taken from Ref. [36] , 



illustrated by Arciprete et al. [134] for the case of GaAs(001)-c(4 x 4). 
The experimental data, which are reproduced in Fig. 14.101 are characterized by 
surface-derived soulder at 0.9 eV (So) and a broad peak at 1.4-2.0 eV (Si), and 
bulk-derived peaks at about 3.5 eV and 5 eV. The latter peaks have also been iden- 
tified in second-derivative off specular geometry spectra at _Eb=100eV by Rowe and 
Ibach [115j as deriving from the bulk critical points, E\ and E^. In the experimental 
spectrum of Fig. 14.101 we have subtracted a background signal, taken to be that of 
the monohydride Si(100)-p(2 x 1):H surface, also reported in Ref. [119J. 

The results of our first principles calculations of HREELS are shown in Fig. 14.101 
for the p{2 x 1), p(2 x 2) and c(4 x 2) reconstructions of Si(100). From the com- 
parison with experiment it is clear that the p(2 x 1) alone cannot reproduce the 
experimental signal since the shoulder at 0.8 eV (So) is missing from its theoretical 
spectrum. A similar observation was made by Gavioli et al. |120j based on tight 
binding calculations of the HREEL spectrum. 

Moreover we observe that, as in the case of the RAS, it is difficult to distinguish 
between the c(4 x 2) and p{2 x 2) calculations. The So peak appears at a slightly 
lower energy (by 0.1 eV) in the p(2 x 2) calculation; however, considering the ap- 
proximations used in the present calculations (scissors shift), it is not sufficient to 
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allow us to prefer the c(4 x 2) over the p(2 x 2). 
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Figure 4.10: REEL spec- 
tra of c(4 x 2), p(2 x 2) and 
p(2 x 1) reconstructions of 
clean Si(100), and compari- 
son with experiment (Farrell 
et al. [119]): E Q = 40 eV; 
6*o = 60°. The surface thick- 
ness is assumed to be d = 
8 layers (plus one vacuum), 
i.e., the half slab. A back- 
ground signal has been sub- 
tracted from the experimental 
spectrum (see text). 
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As discussed in Section 12.2.41 the three-layer model of energy loss is based on 
the assumption that the electron does not penetrate the surface, so that losses occur 
from scattering off long range potentials above the surface. 

In reality, electrons with a 40 eV kinetic energy actually penetrate the surface by 
several atomic layers before elastic scattering occurs. As a result, features in the 
loss which arise from scattering within the crystal itself (as occurs naturally in 
transmission EELS) are missing from our theory. Hence the calculated lineshape 
differs significantly from the experimental one above 2.5 eV. 

To counteract this deficiency of the theory, we augment the reflection loss term 
with a second loss term that represents the trasmission loss, or "bulk" loss within 
the subsurface layers: 

g(q,u) =g 3L (q,u;) + K\q\— (4.5) 
The result is shown in Fig. 14.111 Since the three-layer model already accounts well 
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for losses in the surface layer (the spectra presented here actually do not change 
much if a loss function of the form — l/e s is used), we find we only need to increase 
the proportion of bulklike losses below the surface to improve the agreement with 
experiment. 

In spite of the agreement reached it is clear from the relative intensity of the 
5o and S± peaks that one single reconstruction cannot reproduce the experimental 
spectra. 




1 2 3 4 5 6 

Energy (eV) 

Figure 4.11: REELS spectra of c(4 x 2) calculated using a mixture of long range 
reflection and short range transmittion loss functions, compared with experiment 
(Farrell et al. pH]): E = 40 eV; 9 = 60°. 



4.5.3 Calculated REELS spectra at low energy 

Although many REEL spectra are present in the literature that study the Si(100) 
surface, only a few high resolution spectra that probe the interband transitions of 
Si(100) are available. 

We reproduce in Fig. 14.121 the HREELS data for low-energy incident electrons 
reported by Farrell et al. |119j (Eq = 7 eV; 6q = 60° specular scattering; probably 
at T=300K) and Gavioli et al. [120J (for a slightly different experimental setup: 
E = 6.8 eV; 9 = 62.5°; T=120-500K). 

In both works, two main structures are identified in the spectra below 2 eV. Farrell 
reports an adsorption edge at 0.4 eV with a peak at 0.75 eV (termed So), and a 
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second edge at 1.1 eV with a peak at 1.65 eV (termed Si). We note that the Si peak 
was identified by Farrell as being the peak observed by Ibach and Rowe [1151 1117] 
at 1.7 ±0.5 eV and by Maruno et al. [135] at about 2.0 eV. In the more recent work 
by Gavioli et al, various HREELS data were reported at temperatures ranging 
from 120 to 300 K, and at different analyzer focalizations (see Fig. 14. 12f) . At low 
temperature a shoulder is found at 0.68 eV, as well as two main structures at 0.9 eV 
and 1.15-1.35 eV. It is not clear, however, if the latter feature corresponds to the Si 
structures reported elsewhere, as any peaks appearing at around 1.7 eV in Gavioli's 
data seem to be obscured by the background noise. Note that the main low energy 
peak occuring in the RAS is at 1.6 eV, as seen in Fig. 14.91 which would appear to 
agree with the peak positions of Farrell et al.. In Fig. 14.131 we report the results of 
our ab initio simulation of this HREELS experiment. As noted for the Eq = 40 eV 
data, it is clear that the p(2 x 1) model does not yield the correct lineshape, as the 
So peak is missing. Both c(4 x 2) and p(2 x 2) structures succeed in reproducing the 
double-peaked structure observed in the experiments. In particular, the So peak is 
well reproduced by the c(4 x 2) model, and it is possible that the shoulder observed 
at 0.68 eV points to the coexistence of some p{2 x 2) on the predominantly c(4 x 2) 
surface. The calculated energetic position of the Si peak is in reasonable agreement 
with the data of Farrell, but not so much with that of Gavioli. However, it is worth 



4.5. CALCULATED SPECTRA FOR SI(IOO) 



75 




Figure 4.13: REEL spectra of p(2 x 1), p{2 x 2) and c(4 x 2) reconstructions of clean 
Si(100), for low energy incident beam {Eq = 7 eV; 6q = 60°). Left: Surface thickness 
d = 4 atomic layers; Lorenztian broadening 5l = 0.006 eV; Gaussian broadening of 
5q = 0.3 eV. Right: Surface thickness d = 2 atomic layers. 

to mention that Farrell's data beyond 1.5 eV can be affected to background and 
are not completly relible because some fit of experimental data has been performed 
and in the original paper the dashed line at that energies give us some doubts on 
the exact position of that peak. For this reason we feel more confident in the more 
recent Gavioli's data. 

4.5.4 Analysis of spectra 

In the following two sections we concern ourselves with interpreting the experimental 
peaks observed in the HREEL spectra below 3 eV for the kinematic setup of Farrell 
and Gavioli at Eq « 7 eV. 

Fig. 14. 141 compares the calculated surface dielectric function with the calculated 
HREEL spectra. 

It is clear that there is a one-to-one correspondance between peaks in the energy loss 
and peaks in the imaginary part of the dielectric function. Hence we can analyse £2 
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to characterize peaks in the HREEL spectra. The experimental spectra feature two 
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Figure 4.14: HREEL (along x,y) spectra (top) and surface dielectric functions 
(x, y, z) for E = 7 eV, 9 = 60°. 

main peaks, termed Sq and S± in the literature |119| . The microscopic origins of 
these peaks is now analysed for the c(4 x 2) reconstruction. Figures 14.151 and 14.161 
show the total oscillator strength P#(k), as a function of k, corresponding to an 
energy window of width 25 centred around a chosen peak energy E: 

Jfe(k) = ^|P,, c ,k| 2 (4.6) 

v,c 

for E — 5 < E ck - E vk < E + 5. 

It is clear that Sq arises from transitions located around the T point for both po- 
larizations (see Fig. I4.15P On the contrary Si arises mostly from transitions along 
Y—Y' direction. 

The location of these transitions with respect to the surface band structure 
is shown in Fig. 14.181 Our bandstructure calculation for the clean c(4 x 2) surface 
compares well with that previously published by Fuchs |113j . In particular, a surface 
state is present at about 0.8 eV above the valence band maximum. 

Finally, we looked at |"0n,k| 2 for the valence and conduction band states taking 
part in the strongest transitions. These are plotted in Fig. 14.171 using Xcrysden 
package [136J. We found that S± derives from transitions between surface states: 
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dangling bonds and sp 2 -hydridized p z orbitals, or 7r to it* orbitals as shown on 
Fig. 14.171 To have a confirmation of this statement an analysis of Fig. 14.161 and of 
the band structure (top of Fig. I4.18P is required. In fact, Fig |4.16l shows that the 
most part of transitions contributing to Si comes from k points along Y~Y' path in 
the SBZ and the band structure clarify which bands are involved. In fact, looking 
at the band structure in Fig. 14.181 it is clear that Si peak comes from transitions 
between the flat bands in the Y—Y' direction (blue arrows). 

On the other hand, So is due to transitions between bulk states at the valence 
band maximum (in T) and unoccupied surface states within the fundamental bulk 
bandgap. A representative pair of states involved in this kind of excitation is shown 
in Fig. 14.171 (top). Moreover, with a similar analysis as done for Si, band structure 
and Fig. 14.151 (red arrow) give a confirmation of the fact that Sq is more symetric 
involving transitions around the T point. 




Figure 4.15: Total oscillator strenght Pe(&)(<lE = O.leV) as a function of k 
contributing to So peak. In particular we considered two perpendicular directions 
of polarization belonging the surface plane: X and Y respectively and we analysed 3 
energy windows centred around the position of the peak in the surface epsilon (see 
Fig J4.14jl . In particular 0.32eV for x polarization (left) and 0.285eV (center) and 
0.44eV (right) for y polarization, namely So,a, So,b and So,c respectively. 

4.5.5 Discussion 

We now show how our HREELS calculations can be used to better understand the 
structure of the Si(100) surface as a function of temperature. Si(100) exhibits a p(2 x 
1) LEED pattern above the order-disorder transition temperature of about 200K. 
Our findings (Fig. I4.13P — that the p(2 x 1) surface alone model cannot reproduce 
the experimental HREELS data — are consistent with current understanding of the 
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Figure 4.16: Total oscillator strenght PE(k)(SE = O.leV) as a function of k 
contributing to Si peak. In particular we considered two perpendicular directions 
of polarization belonging the surface plane: X (top) and Y (bottom) respectively and 
we analysed 3 energy windows centred around the position of the peak in the surface 
epsilon (see Fig l4.14p . We considered 3 energy windows centred at 1.1 leV (top), 
for x (left) and y (right) polarization respectively, and at slightly highier energies 
(bottom) 1.18eV (left) and 1.29eV (right) for y polarization, namely Si^S^b, S± t c 
and S^d respectively. 

surface at this temperature (rapid dimer flipping). At lower temperatures, the 
surface should be predominantly c(4 x 2). 

The experimental dependence of the HREEL spectra of Si(100) was studied by 
Gavioli et al. |120j . and their results are reproduced in Fig. 14.121 In that work, 
the observed spectrum was explained as being due to a mixture of symmetric and 
asymmetric p{2 x 1) structures. However, there is much evidence to show that sym- 
metric dimers do not exist at 150K. Furthermore, they did not consider the c(4 x 2) 
structural model, which we find is sufficient to produce a double-peaked HREEL 
spectrum (Fig. 14.131 and I4.10j) . 

Increasing the temperature, the So peak is found to decrease by 5% and the Si peak 
by 40%. As described previously, we took care to ensure that the calculated relative 
intensities of the peaks are not artefacts of our calculation, as we checked their con- 
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Figure 4.17: Isosurface plots of IV'nkMI 2 for represetative states involved in tran- 
sitions responsible for peaks in the low energy HREELS peaks So at 0.8 eV (top) 
and Si at 1.6 eV (bottom). Plots are obtained using Xcrysden [136] package. 
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Figure 4.18: Band structure of the c(4 x 2) (left) with indication of surface peaks. 
Compare with p(2 x 1) (right) to show absence of surface band. A rigid scissor shift 
of 0.5eV has been applied to unoccupied bands. 

sistency with respect to the lorentzian and gaussian broadenings, the surface layer 
thickness, detector integration settings, k-point sampling. We can only conclude, 
therefore, that the c(4 x 2) is not the sole reconstruction present at the surface. 

Furthermore, we note that the p(2 x 1) spectrum does not contribute to the 
Sq peak, on the contrary c(4 x 2) and p(2 x 2) have both some structure at that 
energies. Moreover, looking carefully at Fig. 14.131 it is possible to link the 0.68 eV 
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shoulder in Gavioli's data at T = 150° to the first peak of the p(2 x 2) structure 
and the 0.9 eV peak, present also in the T = 300° curve, to the c(4 x 2) structure. 
At the end we mention that the position of the Si peak needs a correction in or- 
der to predict Gavioli's experimental value between 1.15 eV and 1.35 eV for the 
p{2 x 2) c(4 x 2) and p(2 x 1) structures. Spectra reported in this work, and the 
subsequent analysis, were carried out within the approximation of non-interacting 
particles (RPA), using the DFT-LDA eigenvalues and wavefunctions. A straight- 
forward scheme for incorporating many-body effects is to apply a scissor operator 
to the unoccupied states, following the recipe of Del Sole and Girlanda [37]. In 
this way we compensate for the well known understimation of the DFT-LDA band 
gap, and partially account for self energy and excitonic shifts in energy. A scissor 
shift of +0.5 eV has previously been determined in other works on Si(100) [67] as 
giving the best agreement with the experimental RA spectra. We also confirm this 
result from a fit to the experimental data (see Fig. I4.9|) . Nevertheless, this value 
may not consistently describe the energetic positions of all surface state features, 
which generally undergo many body corrections different from bulk ones. In order 
to determine the correct correspondance between surface-related experimental and 
theoretical energy loss peak, we performed some preliminary calculations including 
many body effects on a smaller (12 layers) c(4 x 2) slab at a lower cutoff (12 Ry). Self 
energy corrections were computed within the so called GW approximation. Within 
this approach it is possible to solve self consistently a closed set of five equations 
(Hedin's equations |137j ) connecting the Green function (G), the polarizability (II), 
the screened Coulomb and vertex interactions (W and T) and the self energy (£). 
with the assumption S = iGW and T = 1. Excitonic and local field effects were 
accounted for by means of solving the Bethe Salpeter equation (BSE), connecting 
the full II with the non interacting II via a 4-points kernel. 

Details of the approach are beyond the scope of this thesis, and can be found in 
Ref. [3]. Our preliminary calculations on HREEL spectra below 2.5 eV show in 
fact that GW+BSE approach give a better agreement with the experimental data 
of Gavioli et al. [120] while the RPA+scissor calculation give a misleadingly good 
comparison with the Si peak of Farrell et al. [119j . Thus the combined GW+BSE 
approach is nowadays the state of the art for computing precise optical spectra, 
nevertheless it is very expensive from the computational point. 
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4.6 Conclusions 

We studied the RA and REEL spectra of the Si(100) surface modeling the surface 
with p(2 x 1) p(2 x 2) and c(4 x 2) reconstructions. Our calculations of RA spectra 
are in agreement with previous works. 

REEL spectra has been calculated for two experimental setup according to the 
available experimental data from Gavioli et al. [120J and Farrell et al. [119J. We 
confirmed that p(2 x 1) cannot be the only reconstruction of the real surface because 
So peak is completely missing in spectra and from a band structure analysis. 
The origin of the So and Si peaks has been carefully analysed considering the c(4 x 
2) model more rapresentative of the surface. We have seen that So arises from 
transitions involving bulk states around T and surface states below the bandgap. 
On the contrary Si involves only surface states. 

Moreover the 0.68 eV feature in Gavioli's data suggest some p(2 x 2) present, along 
with c(4 x 2) but the experimental analysis also suggests that temperature can 
largely change the structural reconstruction because termal motion can easly induce 
a flip-flop of the dimers. 

In summary, we obtained several informations on the nature of the low energy 
excitations of the Si(100) surface and the joint theoretical and experimental REEL 
spectroscopy contributed to clarify the structural composition of this still debated 
surface. 
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Chapter 5 

The oxidized Si(100) surface 



There are two reasons why the investigation of REEL spectra of oxidized Si(100) 
is important. First of all, studying the changes in the experimental spectrum that 
occur after absorption of a foreign substance, such as oxygen or hydrogen, is a 
widely used technique for elucidating the character of spectroscopic features in the 
clean surface. From a theoretical point of view, it is not immediately apparent 
how spectral features related to surface states are modified following atomic scale 
modifications. 

Secondly, a thorough understanding of the oxidation process on Si(100) at the atomic 
scale is of huge technological importance for the development of electronic devices 
and nanodcviccs [138J. 

In spite of an extensive study over several decades, including electron energy 
loss I117j photoemission, RAS \1'69\ EH ED] and several theoretical investiga- 
tions [6? 1 11121 IH3j . there remains some controversy about the reaction pathways, 
with different works suggesting dimer breaking, insertion of O into dimer backbonds, 
as well as silanone bound (0)Si=0 formation |32} I64j. 

In the present section we aim to identify "fingerprints" in the REEL spectrum 
that can help to distinguish between different adsorption sites during the initial 
stages of oxidation, and to obtain further information about the surface states of 
the clean surface. 

5.1 Atomic structure 

Among the vast experimental and theoretical work on oxidized Si(100), a number of 
recent theoretical studies (67J 1112} 1113] proposed various possible adsorption sites at 
low and intermediate coverages, and investigated the optical (RAS) and electronic 
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Structure 


dimer bridge [A] 


dimer backbond [A] 


c(4 x 2) ID 


1.61 


1.61 


c(4 x 2) IE' 


1.78 


1.75 


Structure 


dimer bridge [A] 


Si=0 [A] 


c(4 x 2) ID Si=0 


1.61 , 1.70 


1.53 


p{2 x 2) ID Si=0 


1.63 , 1.73 


1.55 



Table 5.1: Structural lengths of the oxidized Si(100) surface referring to the config- 
urations studied. All calculations are performed in LDA. ID refers to oxygen in the 
bridge, IE' refers to oxygen in the backbond, Si=0 refers to the silanone structure 
created on the top silicon atom of the dimer. 

properties of these systems. 

In this work we considered a wide range of oxygen adsorption sites that we can 




ID IE' Si=0 *o 



Figure 5.1: Schematic diagram of the IE', ID models of the oxidized Si(100) surface. 
Large circles indicate "up" silicon dimer atoms; small filled circles are oxygen atoms. 
The c(4 x 2) unit cells is indicated by shaded regions. 

group according to the number of atomic oxygens involved: those with two atomic 
oxygens in the elementary cell will be referred to as low-coverage (0.5ML) structures; 
on the other hand, a coverage of 1ML will be described in our model by four oxygen 
atoms for cell. 

Among the possible low-coverage structures, we choose to compute the REELS 
for three characteristic local structural motifs (nomenclature follows that of Ref. [1 12|, 
I113j ): structure ID, with an oxygen inserted into the surface Si dimer, and struc- 
ture IE', in which oxygen inserts into the dimer backbond and a silanone structure, 
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whereby oxygen is bonded via a double bond to the silicon dimer. A schematic 
representation of these structures is given in Fig. 15. 11 Both surfaces have a low 
total energy (see Table 15.21 and Refs. |92l 11131 1140] ) and for this reason they are 
good candidates to represent the real system. Regarding the silanone structure, we 
studied a configuration where oxygen is bonded to one silicon atom, as proposed by 
Hemerick et al. in [32] showing a critical formation of silanone species (0)Si=0 dur- 
ing initial oxidation. Anyway, we found silanone structures being higher in energy 
when only two atomic oxygens for cell are considered (see structure "Si=0" with 
just two oxygens added to the clean c(4 x 2) base on top of a silicon atom of the 
dimer in Table 15. 2|) . 

Moreover, in the case of 1ML coverage, we considered several structures appeared 
in litterature, starting from both ID or IE' base configuration. First we studied the 
ID base adding the oxygen on the backbond, resulting a c(4 x 2) 1D+1E' structure. 
Furthemore we considered one oxygen in the bridge and one added in the silanone 
bound because suggested in Refs. [32], [M], finding a relaxed structure with two 
silanone species, one free and the second connected to another silicon by a dative 
bond. We analyzed this "ID Si=0" structure built on both the c(4 x 2) and p(2 x 2) 
bases. The latter case is considered because the two oxygens on silanone bounds of 
adjacent dimers on the c(4 x 2) reconstruction are very close (see black circles in 
the central picture of Fig. 15. 2p . This is unlikely to occur in nature, because there 
would be extra strain on the surface since steric interaction pushes away adjacent 
oxygens. The p(2 x 2) Si=0 structure (see Fig. I5.2p seems to be more reasonable 
for these reasons. 

Table E2] shows a summary of the studied structures with corresponding total ener- 
gies calculated. In order to compare the results, in case of the clean reconstructed 
surfaces we calculated the surface energy: 



where N is the number of atoms in the slab, S is the surface cell area, .Etot is the 
total energy and E'bulk the total energy per atom of bulk Si. The surface with the 
smallest surface energy is the most stable one. 

In case of the oxidized surface we reported the adsorption energy, i.e. the energy 
gain for adsorbing atom at the surface. This quantity is calculated as: 



where Etot is the total energy of the slab, Edean is the total energy of the correspond- 
ing clean surface, Ef ree (0) is the energy of the free oxygen and Nq is the number of 
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oxygen atoms. 

For all structures the inner four layers are fixed to the bulk positions, assuming 
they are not influenced by the surface distortions. Optimization is hence performed 
on the outer slab layers using the Broyden-Fletcher-Goldfarb-Shanno minimization 
(BFGS) |14H 11421 r[l3| 1144] . Crosschecks of results have been performed with PWscf 
and ABINIT codes |123U122] . Because of complexity, the case of one oxygen in the 




c(4x2)-1D-1E' c(4x2)-lD-Si=0 p(2x2)-lD-Si=0 

Figure 5.2: Schematic diagram of the Si-0 related models of the 1ML oxidized 
Si(100) considered in this work. Left: c(4 x 2) with backbond and dative bonded 
silanone; Middle: c(4 x 2) bridge bond and silanone; Right: p(2 x 2) bridge bond 
and silanone. Large circles indicate "up" silicon dimer atoms; small filled circles are 
oxygen atoms. The c(4 x 2) and p(2 x 2) unit cells is indicated by shaded regions. 

backbond and one added in the silanone bond on top of the lower silicon atom of 
the dimer, needs a separate discussion. In spite of the fact that silanone structures 
have been suggested in the litterature (see Refs. [MIES])) both from theoretical and 
experimental grounds, in the present work the total energies founded are relatively 
high for the structures reported up to now. Hence we tried to identify a low en- 
ergy structure for a 1ML coverage that contains the silanone motif. Although there 
are specific techniques for doing such simulations, such as Car-Parrinello molecu- 
lar dynamics, nudged elastic band simulations or potential energy surface mapping, 
we was able to identify one such structure by following a straightforward BFGS 
relaxation, starting from an undimerized Si(100) surface with Si=0 bonds in the 
vertical plane. The eventual purpose is to identify any fingerprints in the REEL 
spectra that might correspond to such Si=0 bonds (see sec. I5.3[) . In Fig. 15.31 the 
evolution of the total energy as this fictitious structure relaxes is shown. We started 
from the undimerized Si(100) with Si-0-Si=0 bonds staying in the vertical plane 
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base 


structure 


E surf [eV/A 2 ] 


N e - 


No 


N fc 


c(4 x 2) 


clean 


0.0903 


256 





6 


p(2 x 2) 


clean 


0.0905 


256 





4 


base 


structure 


Eads [eV] 


N e - 


N 


N fc 


c(4 x 2) 


ID 


-7.14 


280 


2 


6 


c(4 x 2) 


IE' 


-7.22 


280 


2 


8 


c(4 x 2) 


Si=0 


-5.73 


280 


2 


2 


c(4 x 2) 


1D+1E' 


-7.17 


304 


4 


2 


c(4 x 2) 


ID Si=0 


-6.76 


304 


4 


2 


p(2 x 2) 


ID Si=0 


-7.04 


304 


4 


4 


c(4 x 2) 


"A" in Fig. E2 


-6.74 


304 


4 


2 


c(4 x 2) 


"B" in Fig. E3] 


-6.82 


304 


4 


2 


c(4 x 2) 


"C" in Fig. 1531 


-7.45 


304 


4 


2 



Table 5.2: A summary of energetics of all stable and metastable configurations 
found after BFGS relaxation for the oxidized Si(100) surface. Structures are listed 
with their total energies, each force component is relaxed below the threshold of 
0.13 eV/A. 

(perpendicular to the surface). After some BFGS steps, the Si-Si dimers are formed, 
but the Si-0-Si=0 bonds are still both present. In configuration "A" (see Fig. 15. 3p . 
one of the Si=0 starts to bond the dimer until the metastable structure "B" where 
there is still the Si-0-Si=0 bond in the plane (see blue arrows in Fig. 15. 3p . but the 
other oxygen is bonded in a Si-O-Si-O-Si chain, similar to a part of the Si02 crystal 
lattice. Leaving this structure to relax we found the configuration "C" without any 
silanone bonds and characterized by a silicon atom bonded to 3 oxygens. In this 
last configuration Si(100):O (1ML) reaches the lowest total energy. 
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Figure 5.3: Total energy of configurations generated by relaxing Si(100) with one 
oxygen on the backbond and one oxygen on top of the higer silicon atom of the 
dimer. 

5.2 Electronic structure and the effect of oxidation 

We present here the computed Kohn-Sham band structure of three models of ox- 
idized Si(100). The electronic structure of the ID and IE' surfaces is shown in 
Fig. l5.4l along a standard path in the surface brillouin zone: T-X-Y'-Y-T (see Fig. 12.21 
in chapter [2]). Mostly, modifications respect to the clean Si(100)-c(4x2) concern a 
change of surface states inside the bandgap. In fact, comparing Fig. !4.18l to Fig. 15.41 
it is evident that oxygens move the surface states visible in the T-X direction, states 
that are responsible of the So EEL peak of the clean surface, (see Fig. 14.131 and 
Fig. I4.12l in Chapter 2]). The oxygen adsorbed modifies also bands in the Y'-Y path. 
Moreover, we present in Fig. 15.51 the computation of the electronic states for the 
surface including a silanone (0)Si=0 bound. In Fig. 15.51 it is possible to see the flat 
bands in the Y'-Y path due to the presence of oxygen doubly bonded to the top 
silicon of the dimer. Furthermore we performed a test on the p{2 x 1) structure to 
better analyse the effect of the isolated Si=0 bond. For this system we report the 
bands compared to that of the clean p(2 x 1) surface (see Fig. 15. 6p . The bandgap 
typical of a semiconductor is closed by the presence of oxygen, hence the system is 
metallic. 

All calculations are performed within the DFT-LDA and a scissor shift of +0.5 eV 
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Band structure Si(1 00) :0 1 Ep 




r XV Y r 

Band structure Si(1 00) :0 1D 




r xv y r 



Figure 5.4: Band structures of the IE' (top), ID (bottom), models of the oxidized 
Si(100) surface. A rigid upward shift +0.5eV (scissor shift) has been applied to the 
unoccupied bands. 



(see sec. I4.5.5P has been applied to the unoccupied bands in order to mimic the 
many body effects. Band structure calculations are performed with the ABINIT 
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Band structure Si{100):O 1D SiO 




r XV Y r 

Figure 5.5: Bandstructure of the ID Si=0 model of the oxidized Si(100) surface. 
A rigid upward shift +0.5eV has been applied to the unoccupied bands. Flat bands 
along the Y'-Y direction are due to the presence oxygen double bonded to the top 
silicon in the dimer. 



package [122] . 



5.3. EXPERIMENTAL ENERGY LOSS DATA 



91 




Figure 5.6: Bandstructure of the clean Si(100)-p(2xl) (black line) compared with 
the ones of p(2 x 1) Si=0 (atomic oxygen) (red dashed line). 

5.3 Experimental energy loss data 

Experimental EEL spectra of the oxidized Si(100) surface are not so numerous in 
litterature, at least in the range of energy we are studying (i.e in general up to ~8- 
10 eV and in particular up to ~3 eV). In this work we are refering to not-so-recent 
data from Ibach et al. I117j and at the subsequent paper of Ludeke et al. |118j 
where experimental REELS are reported in terms of the second derivative of the 
spectrum respect to the energy in order to cut the large contribution of the elastic 
peak and evidence spectral features. 

In Fig. 15.71 all experimental data are collected. The spectra from Ibach et al. are 
reported for the clean (black circles) and oxidized (red circles) surface. In particular, 
Ibach used high incident energy of the electrons, (Ej=100 eV) and an off specular 
geometry with normal incidence and collecting electrons at 42.5° respect to the 
perpendicular direction to the surface plane. 

From the analysis of Ibach's data we can conclude that there are two features 
that seems to be due to the oxygen adsorption: a peak at around 7 eV, which 
appears in the oxidized surface, and the surface plasmon at around 12 eV which 
is splitted in two. Moreover, the bulk plasmon, appearing at 17 eV in the clean 
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Figure 5.7: Experimental data from Ibach et al. and Ludeke et al. for the clean 
(black) and oxidized (red and blue) Si(100) surface. 



surface, is lowered and slightly blueshifted. 

The other experimental data reproduced in Figure 15.71 are from Ludeke et al. 
|118j (green, blue and violet circles in Fig. I5.7P and refer to specular scattering 
geometry: backscattered electrons are detected and analyzed. The three sets refers 
to different incident energies: 24 eV (violet), 50 eV (blue) and 100 eV (green). Peaks 
at 3 eV and 5 eV (Ei and E2 respectively) appear in all those sets of data. From 
the comparison of the three sets it is possible to conclude that the two main peaks 
around 5 eV and 7 eV are in agreement with Ibach's results. Moreover the 7 eV peak 
increases with increasing incident electrons energy and is interpreted by Ludeke et 
al. as an excitation of the SiO or Si02 molecules |145l 1146] in the Schumann region. 
We will comment this assignement in the following paragraphs, where we will show 
that our results do not support it. 



5.4 Theoretical energy loss spectra 



In the following paragraphs we present REEL spectra calculations of the oxidized 
Si(100) surface using the code YAMBO [126] . 
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5.4.1 Clean versus oxidized Si(100): the low-energy part of the 
spectrum 

REEL spectra are often presented in literature in terms of the second derivative 
of the energy loss. Moreover we start presenting here our computed bare REEL 
spectra for the clean and oxidized Si(100). In particular we considered the clean 
and the 0.5 ML oxygen-covered surfaces with oxygen inserted into the bridge (ID) 
or the backbond (IE'). We compare the REEL spectra computed for the experi- 
mental setup of Farrell (Ref. [119j : specular scattering at 60° with incident energy 
Ej=40 eV) with measured spectra for the clean Si(100). All results are shown in 
Fig. 15.81 where theoretical spectra for x and y polarization directions are also pre- 
sented separately. In the bottom panel of Fig. 15. 81 a comparison of the (unpolarized) 
experimental data with computed spectra averaged over the two polarization direc- 
tions is shown. 

It is evident that the two main peaks, namely the So shoulder at 0.8 eV and the peak 
Si at 1.3 eV in the experiments, appear also in the calculated spectra for both the 
clean (black line) and the oxidized IE' configurations (blue line). On the contrary, 
oxygen in the bridge position (ID structure-red line) completely kills the first peak 
of the spectrum. Unfortunatly we did not find experimental data for this surface at 
low oxygen coverage in the low energy range (< 3 eV), but a similar behaviour is 
shown in the case of a coverage of H2O in Ref. |119j . In conclusion of this section 
we can say that REEL provide a useful tool able to distinghuish between oxidized 
surface configurations, in particular in the low energy spectral region. 

5.4.2 Clean versus oxidized Si(100): REEL spectra in a wider spec- 
tral reagion 

In this section we present computed REEL in terms of the second derivative spec- 
trum. In Fig. 15.91 (top panel) we show the bare REEL spectra for the clean, the ID 
and IE' structures (0.5 ML coverage) computed for the experimental setup used by 
Ibach and Ludeke in their works |116 I118j (i.e. incident energy Ej=100 eV with 
slightly different geometry scattering^]). In the middle panel of Fig. 15.91 we show 
the computed second derivative spectrum. Neglecting oscillations for energies below 
3 eV, and except for the intensity of the surface plasmon peak calculated around 



1 Ibach used off specular experimental geometry instead of backscattering geometry form a normal 
incident beam of electrons used by Ludeke el al. in Ref. 118 . However, at this level of approximation 
we did not find appreciable differences in case of spectra calculated with specular or off specular 
geometry. 
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Figure 5.8: Experimental and computed REELS at Eq = 40 eV for clean Si(100)- 
c(4 x 2) is compared with the calculated spectra of ID, IE' and ID Si=0 oxidized 
reconstructions. Experiments for the clean are taken from Farrell et al. in Ref. [119J. 



15 eV, the clean and the other two oxidized surfaces do not show any important 
differences. This is in contrast with experiments (bottom panel of Fig. 15. 9p where 
a peak centred at around 7 eV appears after oxidation. In order to understand the 
origin of this peak, we analysed the EEL spectra of all configurations presented in 
Section 15.11 Following arguments presented in Ref. |118| about the origin of the 
7 eV peak, we show theoretical results for surfaces including a silanone bond on the 
c(4 x 2) and p(2 x 2) reconstructions of the clean Si(100) and an oxygen inserted 
in the dimer bridge (Fig. 15. lOj) . We found slight differences with respect to the 
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Figure 5.9: Computed bare 
REEL spectrum (top panel) 
for clean Si(100)-c(4 x 2) com- 
pared with the ID (oxygen in 
the bridge) and IE' (oxygen in 
the backbond) oxidized recon- 
structions. Calculated neg- 
ative second derivative spec- 
tra of the same surfaces (mid- 
dle panel) shows slight dif- 
ferences in the spectral re- 
gion above 3 eV. Experimen- 
tal data (bottom panel) are 
also reported from Refs. |117] 
and [118j showing an impor- 
tant peak at around 7 eV, be- 
yond the main peaks at 3 eV 
and 5 eV, interpreted as the 
Ei and E2 bulk silicon transi- 
tions. The 7 eV peak is absent 
in the experimental spectrum 
of the clean (black circles) and 
is hence interpreted as due to 
oxidation (red and green cir- 
cles). 



previous calculations: the 7 eV peak still does not appear and the spectral features 
are similar for the clean and the oxidized surfaces. Unfortunatly, large differences 
between c(4 x 2) and p(2 x 2) structures appear in the low region of the spectrum 
where experimental data are not available. However we can observe that the surface 
plasmon is slightly lowered and blueshifted when a silanone bound is included in the 
surface configuration. Moreover we report in Fig. I5.10I (left) the results concerning 
the stable "C" structure and his metastable precursor "B" containing one (0)Si=0 
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Figure 5.10: Left: REELS for clean Si(100)-c(4 x 2) compared with the ID Si=0 and 
the p(2 x 2) ID Si=0 reconstructions. Right: REELS for clean Si(100)-c(4 x 2) com- 
pared with the B structure (metastable precursor) and C structure For comparison, 
experimental data from Refs. |118t I117j are reported. 

bound depicted in Fig. [57 



Once more we can see that no features at 7 eV comparable to the experimental 
peak appears in our calculations, and even if Hemerick et al. [32] suggests similar 
structures on small clusters as energetically favorable, we can not conclude that 
silanone bond or siloxane structures on Si(100) are representative of the real sur- 
face. Moreover, Ludeke et al. guessed that the 7 eV peak was related to molecular 
excitations of Si=0 bond or silicon monoxyde molecule. Calculations performed in 
the present work, however, rule out Si=0 structures bonded to the surface. 



2 For c(4 x 2) silanone structures we used 72k points in the IBZ and 450 bands, corresponding 
to minimum e—h transition energy of 13.27 eV p(2 x 2) required 64k points and 500 bands, corre- 
sponding to 15.09 eV, for the "B" and "C" structures 72k points and 500 bands, corresponding to 
15.09 eV, have been considered in order to fully converge. 
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It can be worth to mention that, taking as a reference the surface plasmon peak 
of the clean Si(100) (represented by the black continuous line in Figs. 15.91 I5.1(jp . 
oxygen added into ID and IE' positions lowers the peak (see Fig. 15. 9p . while oxygen 
added into the silanone bonds lowers and redshifts the position (see Fig. 15.101 left 
panel). Conversely, calculated spectra for B or C configuration show a shift of the 
peak position to the opposite direction (see Fig. 15. 1UI right). 

At this point, several speculations can be done in order to explain the origin of 
the peak measured at 7.0 eV, which does not appear in theoretical spectra for the 
considered structures. From the point of view of calculations we must underline that 
they are performed only at the RPA level. A further analysis could call for the use of 
approximations beyond RPA in order to describe eventually strong excitonic effects. 
Moreovoer, we mention that the oxidation process of Si(100) is still under debate, 
for example, we can not exclude the formation of clusters of SiC>2 during oxidation. 
In addition, looking at the intensity dependence of the peak, increasing with the 
energy of the incident electrons (see Fig. 15 . T[) . we could interpret the peak as due to 
transitions from states originating from structures below the surface level, in fact, 
higher energy electrons are expected to penetrate deeper in the sample. Within 
these hypotesis the present theory could not be adeguate to treat, for example, 
multiscattering processes. In addition we note that the experimental data are old 
and it is not clear how well characterized and clean the surfaces are. However we 
conclude that our calculations does not support Ludeke's interpretation about the 
origin of the 7 eV peak, because does not appear in all the considered model surfaces 
(including Si=0 or O-Si-0 bonds). 



5.5 Conclusions 

In the present section we draw a summary of the previous analysis, and some conclu- 
sions. We calculated the relaxed atomic position of several oxidized si(100) surface 
reconstructions at 0.5ML (ID and IE') of coverage. Moreover we analyzed 1ML 
coverage structures with a silanone bund and an oxygen in the bridge of the dimer 
(ID Si=0), with a p(2 x 2) and a c(4 x 2) base. Finally we found a structure "B" 
characterized by a silanone bond and Si-O-Si-O-Si chains, to be the metastable pre- 
cursore of the "C" structure, the most stable, with one silicon atom bonded to three 
oxygens (see Fig. 15. 3|) . 

We calculated the bare REEL spectrum in the cases of the oxydized Si(100) pre- 
viously considered. We have been able to relate the change in the spectra with the 
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changes in the band structure in the case of both prototypical configuration (oxygen 
in the bridge and oxygen in the backbond). Furthermore we calculated the second 
derivative spectra thanks to the implementation of the SG method in the YAMBO 
code [126] . This allowed us to compare EELS results with the experimental data. 
We have shown that the excitation at 7 eV can not be related to a molecular exci- 
tation of Si=0 or O-Si-O, as suggested by Ludeke et al, because the configuration 
with an oxygen on top of the higher silicon atom in the dimer does not reproduce 
that structure. At the moment we are not able to give an alternative explanation 
for the origin of this spectral feature, but we can confirm that the reconstruction 
studied in this paper (including the most stable structures ID, IE', ID Si=0 p(2 x 2) 
and c(4 x 2) "C" and the metastable precursor with Si-0 and Si=0 molecule) do 
not reproduce this peak. Still, recent works by Hemerick et al. and Chabal et al. 
have pointed out the critic presence of silanone bonds during the oxidation process. 
We can only guess that including many body effects in the calculations could pro- 
vide more accurate theoretical results, taking into account self-energy and excitonic 
effects, in order to help in shedding light on that problem. 



Chapter 6 



Subtleties in electronic 
excitations of open shell 
molecules 

In this chapter we present a theoretical study of BeH, a simple heteronuclear di- 
atomic molecule with an unpaired electron. We considered BeH molecule as the 
simplest exemple of an isolated open shell system for which the ground state is ex- 
pected to be spin polarized. We present calculations of energy levels and density of 
states. Convergence studies and problems arising up in order to correctly describe 
excitation spectra are also underlined. At the end some excitation energies calcu- 
lated for this system are compared to available experimental data and the agreement 
is discussed. 

6.1 Brief review of TDDFT for isolated systems 

In this section we briefly review the Casida's approach to treat the molecular exci- 
tation of isolated systems. This framework is a reformulation of the TDDFT in the 
configuration space and the details can be found in Ref. [147} 1148] . In fact all the 
observables are represented using the DFT-KS eigenvectors in such a way that the 
non-interacting response function results to be diagonal. 

We note that in case of isolated systems the physical quantity we deal with is the 
polarizability of the system a(uj) that is defined when an external perturbation 
8V ext (t) = z5E z (t) is applied to the system. In this case the x component of the 
dipole momentum is defined by: 5d x = —q5x where q is the charge of the electron 
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and 6x = J d 3 x5p(x,t)x. 

Hence, the polarizability is defined by: 



d 3 x 



E z (x) 



and using the Lehmann representation it can be written by: 

(Wo|£|*j}(¥j|2|*o> 

i 



a. 



( El -Eo)*-u,* 



•1) 



(6.2) 



Where Ej, Eq, \&o are the energies and the many body wavefunctions of the 
excited and ground state respectively. Using the basis set of operators {alu? ,a ma } 
we can define the quantity: 



P ija (t) = (* KS \al(t)a ja (t)\V 



,KS\ 



(6.3) 



corresponding to the density in the configuration space, for which the following 
relations held: 



SP ij a(h) = f dt>x?£ hkT (h-t 2 )5V hkT (t2 



5Pija(tl) 
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(6.4) 
(6.5) 



where the KS response function is written over this basis set: 
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In case of independent particles and applying the a and at operators we obtain the 
reduced expression for the response function: 

fja fia 
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Hence, in the configuration space the TDDFT equation held the form (see Ref. [147 
EH]): 



klr 



x x x u ~ ( gfeT ~ £kr ^ IT (, A 

OarOjhOik 7 ~? J^ija,hkr\ u ) 

Jkr JHt 



5P hkT {u) = 5V*?{u) (6. 



This equation can be reorganized and written as a matrix equation (see Ref. [1471 
1148] for mathematical details) exploiting the hermitian nature of the kernel Kij^kr- 
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Moreover assuming that the KS orbitals and the kernel are real quantities (infinite 
lifetime of the excitations) and after some algebra we can write: 



A(u) + B(u) 
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and the polar izability as: 



OarOihOjk ~? _ f 

Jhk Jkr 

—Kijcr^hri 00 ) 
fhr fkr 



Kija^hkriu) 



(6.9) 
(6.10) 
(6.11) 

(6.12) 



where S(u) = C(A-B)~ 1 C and O(w) = S- 1 / 2 {A+B)S- 1 / 2 . Finally, comparing this 
last expression (|6,12p with eq. (|6.2p it is possible to write the eigenvalue equation: 



where: 

^(wz)ijV,ft*r( w ) 



Si,hSj,kSa,r(£k,cr ~ ^huf + 



.13) 



(6.14) 
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Eq. (|6.13|) is called Casida equation (see Ref. [147} I148| ) and allows to map the 
problem of finding the excitation spectrum of a system too an eigenvalue problem 
for the matrix Q(ui). 



6.2 Closed and open shell systems: the problem of the 
correct counting of excitations 

We divide now the discussion to the case of closed and open shell systems. In the 
first case the approach presented in the previous section predicts, in principle, the 
correct excitation energies. However, in the second case some problem limits its 
straightforward applicability. In the next paragraphs we present the problem with 
an example. 
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6.2.1 Closed shell systems 

Let us start with the case of the closed shell systems, i.e. systems having the same 
wavefunction for the two spin channels. In this case it is possible to diagonalize the 
matrix, presented in the previous section, exploiting the symmetries resulting from 
the adiabatic approximation. Let us consider a three levels system consisting of two 




-H- 



- 




4- 



Figure 6.1: Schematic representation of the possible excitations generated in a 3 
levels system with 2 electrons with opposite spins. 

electrons (see Fig. I6.1|) . the ground state is characterized by two coupled electrons 
both occuping the i level, two additional unoccupied levels j and k are available to 
the excitations. A spin-unpolarized ground state is expected and the £1 matrix has 
important symmetry properties, so that its eigenvectors will be either symmetric 
or anti-symmetric with rispect to spin flip. Assuming the adiabatic approximation, 
the f2 matrix, can be written as: 

such that : 



fyt = and n n = %t 
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The eq. f)6.13f) is decoupled in: 

F^ = F n , (n n + n n )F n = ujjF It (6.15) 
F n = -F n , (n n -n n )F n = ujF n (6.16) 

In case of eq. (|6.15p for a given electron-hole pair, both spin channels interfere con- 
structively, TD total charge changes, but not the TD magnetization (Spin-singlet 
excited state). Conversely in case of eq. (|6,16p both spin channels interfere destruc- 
tively, no TD total charge change, TD magnetization change (Spin-triplet excited 
state). In this case the resulting eigenvectors are written in the form: 





Hence, summarizing, it is possible to group the results in two kind of excitations: 



1 

71 



4>l /2 = 4= (|j t U) + H 1 3 I)) singlet (6.17) 



and 



I* t 3 I) (6-18) 
4 = -L(|jtU>-|itjl» triplet (6.19) 

<t>\ = (6.20) 

This states are simultaneous eigenvectors of H KS , S 2 and S z . 
6.2.2 The problem of the open shell systems 

The simplest case of a system with spin-polarized ground state consists of three 
electrons that can be arranged within three levels. Fig. (|6.2p describe the ground 
state configuration and the possible spin-collinear excitations. In this case the Q 
matrix does not have the same symmetries presented in the previous section. This 
time it contains the elements related to the excitations jk f, ij \,, ik f and ik \, 
resulting with 4x4 dimension. Moreover any symmetry does not allow to divide 
by blocks the whole fi. However the selection rules allow to group the excitations 
in two sets: conserving the spin AS" = and those for which AS = 1. 

Diagonalizing H KS , S 2 and S^ simultaneously we obtain the following doublet 
for the ground state: 
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Figure 6.2: Schematic representation of the possible excitations generated in a 
3 levels system consisting of 3 electrons. The yellow square represents the double 
excitation missed by the casida's framework. 



and the excited states (doublet): 



and the quadruplet: 



&1/2 
*l/2 

% 



\iik) (6.21) 

(6.22) 

i= [\ijk) + \ijk)) (6.23) 
i= + - 2|tjfc)) (6.24) 



+ \ijk) + (6.25) 



where we adopted the notation i to indicate the \, spin occupation of the i state. The 
double excitation \ijk), represented by the yellow square in fig. 16.21 is present because 
the eigenvectors must be simultaneous solutions of the eigenvalue problem for H KS 
and S 2 . The excitations of this system present states with three half occupied levels 
(see Fig. I6.2p that are mixed by the S 2 operator with similar states with inverted 
spins. 
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It is worth to note that this kind of excitations are missing when we attempt to 
build the eigenstates of H KS , S 2 and S 2 starting from the single particle excitations. 
In this case we obtain three doublet: 

(6.26) 
(6.27) 

k) + \ijk)) (6.28) 
and a triplet: 

4>\ /2 = -j= {\ijk) + \ijk)) (6.29) 

As the theory is exact and the sole quantity requiring an approximation is the 
kernel the lack of the excitations description is due to the adiabatic assumption. 
The research of an improved kernel able to correctly describe the missed excitations 
is beyond the scope of this thesis. 

However, in the following paragraphs an application to the case of the calculations 
of the excitation energies of the BeH molecule is presented. 



6.3 A simple open shell molecule: BeH 

A BeH molecule consists of five electrons, four due to Beryllium and one to Hy- 
drogen. Configuration of the Beryllium atom ( 4 Be: Is 2 2s 2 ) suggests his divalent 
nature, in fact usually the two electrons in the Is level are not involved in chemical 
bonds. Therefore we describe the joint effect of the core and the inner electrons 
of Beryllium with a suitable pseudopotential. In particular we used a TrouUier- 
Martins LDA pseudopotential available in the on line repository of the ABINIT 
code [681 1122| . Following this recipe, only three valence electrons of the molecule 
are involved in DFT calculations and are expected to determine the main properties 
of BeH, including the spin polarization of its ground state. 

Bond length between Beryllium and Hydrogen is assumed to be R e(? = 1.3426 A 
(see Ref. j!49 j) and the cubic supercell is centred in the middle point of the line 
joining the centers of the two atoms. A schematic representation of the molecule is 
depicted in Figure [6731 where the active electrons (blue and red circles) are underlined 
and distinguished respect to the spectator electrons (white circles). Spin variable 
is also shown in order to evidence that spin is not compensated. In the following 
paragraphs we will show our results about spin resolved calculations performed with 
the ABINIT package [68] . 
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Figure 6.3: Schematic view of the atoms in the cell used to represent the BeH 
molecule (acell is not on scale). Red and blue circles represent active electrons, 
white circles represent Beryllium Is electrons. The effect of these core electrons is 
described by a suitable pseudopotential (see text). 

6.4 Energy levels and density of states 

In order to obtain converged energy levels and the spin resolved density of states 
of BeH, we performed total energy minimization within DFT in local spin density 
approximation. 

We obtained the theoretical energy levels (see Fig. 16. 4p imposing a conver- 
gence tolerance of 10 15 Ha 2 for the largest squared residual defined by: (nk\(H — 
E) 2 \nk), where E = (nk\H\nk) . In order to establish reliable results we performed 
convergence tests on the supercell size (see Sec. 16. 4. Tj) . Here we comment the results 
summarized in Table 16.11 We compared our results with theoretical calculations 
reported in Ref. |149| and experiments from Ref. [150]. We obtained a satisfactory 
agreement of the Kohn and Sham energy of the highest occupied orbital with respect 
to calculations performed in Ref. |149] with localized basis orbital. 

Even if there is not analogous of the Koopmans theorem, providing a physical 
interpretation of the Hartree Fock eigenenergies, for the KS eigenvalues, it is possi- 
ble to interpret the highest occupied molecular orbital (HOMO) energy calculated 
within DFT as the ionization energy of the system (see ref. [151] ). 
However in the present case of the BeH molecule, if we assume that the E3 CT . energy 
referred to the vacuum should be the ionization energy of the system, we face with 
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Figure 6.4: Spin resolved Kohn-Sham energy levels for the BeH molecule calculated 
within DFT-LSDA. Energy levels are calculated for a 55 Bohr cubic unit supercell 
with a plane-wave cutoff energy E cu j=18 Ha and 20 Kohn-Sham states.. 
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Figure 6.5: Spin resolved density of states for BeH calculated within DFT-LSDA 
with the code ABINIT [BH]. Red and blue lines represent the spin up and down 
components. A smearing of 0.008 Ha is applied and a 55 Bohr length unit cell is 
considered. 

a discrepancy with respect to the experimental value reported in Ref. [150J. In fact, 
the experimental value of 8.21 eV disagrees with - £ homo _ 462 eV) calculated 
within DFT-LSDA, and also with 4.60 eV obtained in Ref. |149j . This discrep- 
ancy is due to the inadequacy of the LDA exchange and correlation potential (see 
Ref. |152} 1153] ) that does not have an asymptotic — 1/r required behaviour. This 
fact lead to the lack in the description of the levels staying close to the vacuum 
region and the excitation energies of the molecule involving these levels. 

Figure [631 represents the distribution of the density of states obtained after con- 
vergence tests, discussed in detail in Sec. 16.4.11 and, in Figure 16. 5( we report the 
density of states (DOS) calculated within the same approximations. When spin 
symmetry is broken, BeH molecule shows an asymmetric distribution of DOS (see 
Fig. 16. 5p with respect to the two spin channels. We considered a cutoff energy 
Ecutoff = 18 Ha and a cell size ao = 55 Ha as converged values. The first three 
peaks, centred at —8.35 eV, —8.01 eV and —4.62 eV represent the occupied states, 
also depicted in Fig. \6A\ moreover 3a^ and In orbitals are shown, in particular 
the height of peaks related to Itt orbitals is doubled with respect to the a orbitals 
bacause of their degenerancy giving contributions to the DOS peak coming from a 
double number of states. A smearing of 0.008 Ha is applied to both distributions. 
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Figure 6.6: Top view of the BeH molecule (black circles) inside a cubic cell, grey 
zone represents the vacuum region around the molecule (cell is not on scale). As 
an exemple we report a comparison between two possible choices of the molecule 
placement inside the cell: distance between atoms of the replicas is larger in case 
(a) than (b). 



6.4.1 Convergence issues 

Convergence tests have been performed respect to two quantities: (i) the cutoff 
energy E cu t determining the number of plane wave of the basis set and (ii) the 
supercell dimension acell. 

In order to determine a suitable value for E cn t we checked the values of the 
KS-eigenvalues: E2cr t , E2o-^ and E3 CTt , corresponding to the occupied states (see 
Figure E3 and Table IBTTj) . Since we used a Troullier-Martin (TM) pseudopotential 
for Be, we found a resonable converged value of E cu t = 18Ha (corresponding to 
36Ry) see Table 16.11 and Fig. 16.41 However we verified that E cu t = 10 Ha can 
be considered a converged value in order to calculate the first excited states of the 
molecule, in fact, the HOMO-LUMO distance for each spin channel does not change 
significantly. Table 16.11 summarize the converged values for A-j-j- = lir-f — 3cr^ and 
A44 = 3(7^ — 2cr^, that distances are also evidenced in Fig. 16.41 Moreover, a set 
of tests with respect to the cell size have been performed in order to find the best 
compromise between accuracy and practical feasibility of calculations. In fact, a 
large vacuum region around the molecule is required in order to obtain the correct 
behaviour of continuum states above the vacuum level (see Fig. 16. 7p and to avoid 
spurious interactions between replicas, however this request can increase significantly 
the computational workload. 

It is worth to mention that a good choice of the cell shape (see Fig. I6.6P can slightly 
optimize the numerical convergence. We considered the molecule placed in the 
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Level 


present work 


Ref. [154J 


Ref. [149] Exp. Ref. [150] 


2cr t 


-8.35 


-8.35 






-8.01 


-7.99 




3cr^ 


-4.62 


-4.63 


-4.60 8.21 


HOMOLUMO 


present work 


Ref [ES] 






2.28 


2.29 






5.07 


5.16 





Table 6.1: Theoretical spin resolved energy levels and HOMO-LUMO distance for 
BeH molecule calculated within DFT-LDA using the code ABINIT [68]. A 55 Bohr 
cell, E cut =18 Ha and 20 bands have been used and considered converged values. 

middle of a cube of length ao (acell) with the straight line joining Beryllium and 
Hydrogen atoms parallel to one face of the cube (Fig. 16.6( b)). 

We performed several tests increasing ao for a fixed cutoff energy and number 
of bands. In Fig. 16.71 a summary of the results is depicted: in the actual range 
of ao considered, occupied states slightly change because of the dimension of the 
supercell, on the contrary the unoccupied levels are strongly influenced and require 
a large cell size (ao > 50 Bohr). This fact can be important if we want to calcuate 
excited energies involving higher energy levels. In Fig. l6.4l we can distinguish discrete 
unoccupied states (l7r, 3c, 4a) and a thick serie of states above the vacuum level (grey 
region). In fact, increasing the cell size, the upper levels become closer and closer 
with each other (see Fig. I6.7P and a continuum of states, in the limit case of an 
infinite supercell. In conclusion we assumed that a cell size of 55 Bohr is required 
in order to correctly calculate the three lowest eigenvalues. On the contrary, for 
excited states calculations, involving higher energy levels in the vacuum, a larger 
cell is recommended. Our conclusions are in agreement with a slightly different 
treatment discussed in Ref. |154| . 
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Figure 6.7: Energies of the molecular orbitals of BeH as a function of the cell size. 
Spin up (right panel, red circles) and down (left panel, blue circle) are distinguished. 
Unoccupied states are strongly influenced to the cell size and require ao > 50 Bohr. 
Tests are performed with a cutoff energy of 10 Ha. 



6.5 Excitation energies 

In this section we report results concerning excited state calculations within the 
TDLSDA formalism. For that we used a recent implementation of the method in 
the ABINIT software [681 Q221 EM]. In Table O we summarize the results about 
the main excitation energies calculated. We did not consider excitations involving 
continuum states because they are extremely influenced by the countour conditions. 
However, in order to have more reliable results we used a larger supercell with respect 
to the case of ground state calculations, in particular a 70x70x70 Bohr cubic cell 
is considered with a cutoff energy of 10 Ha and 50 bands. 

Even if the ionization energy disagrees with respect to experimental data making 
a shorter distance between occupied states and vacuum level, our first calculated 
excitation energy n connecting the spin-collinear states 3a —> In is comparable with 
experimental value reported in Ref. [155] . The TDLSDA theoretical value (2.39 eV) 
understimates the experimental measure (2.48 eV) but improves the DFT-LSDA 
HOMO-LUMO distance (2.27 eV) calculated as a bare difference between Kohn- 
Sham eigenvalues (see Table 16.11 and I6.2[) . 

However we observed that every excitation involving states higher than the 4cr 
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Excitation 


elec. conf. 


this work 
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Table 6.2: Theoretical excitation energies for BeH molecule calculated within the 
TDDFT framework and ALDA approximation using the code ABINIT [BH] . Energy 
values are reported in eV. Converged results are obtained with a cutoff energy of 
E cu t=10 Ha, 70x70x70 Bohr and 50 bands. We compare our results with calculated 
values within two TDDFT implementations in the code DeMon2K (see Ref. [149, 
I156j ). Exp. a are reproduced from Tab. I of Ref. [157] and Exp. b from Ref. |155j . 



level are strongly influence by the cell size. For this reason we evidenced with bold 
carachter in Table 16.21 the more reliable results, i.e. excitations excluding states 
above 4<r. Moreover we underline that the II excitations, involving three half occu- 
pied states (2o" 1 3(T 1 l7r 1 ), appear at two excitation energies (5.66 eV and 7.15 eV) and 
with the same orbital configuration. These Il-excitations are distinguished by their 
total spin, i.e. S= 1/2 in one case and higer total spin in the other. Thanks to the 
Hund's rule we can identify the excitation with S= 1/2 with higer energy (7.15 eV), 
and this lead us to discard the 5.66 eV excitation. A different argument leads Casida 
et al. to reject the same excitation identifying the same excitation we discarded as 
spin contaminated. Comparison with experiments available in Ref. [155| I157] are sat- 
isfactory, in particular with respect to the case of II : 3cr-|- — > l7i> and II : 2<r^ — > l7r+i 
excitations. Discrepancy about II : 3cr^ — > 2tt^ excitation is due to the problem pre- 
viously underlined: in fact 2ir states are higer than 4a level. 

Comparison with results obtained in Ref. (14=91 1156] are also aligned with our con- 
clusions. 
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6.6 Conclusions 

We studied the BeH molecule as the simplest exemple of open shell system. In our 
pseudopotential approach only three valence electrons are involved in the eigenvalue 
problem, for this reason total spin is not compensate and the electronic ground state 
is spin polarized. 

Spin resolved density of states and energy levels have been computed within 
DFT-LDA approach. Convergence tests with respect to the supercell size have 
shown the sensitivity of the unoccupied states above the vacuum level to the cell 
dimension. Moreover we found a discrepancy between the energy of the HOMO or- 
bital calcualted (in agreement with others theoretical works, see ref. |149] and [156 ) 
and the ionization energy measured. This is due to the well known inadequacy of 
the asimptotic behaviour of the LDA exchange and correlations potential. 

In conclusion we calculated the first excitation energies of the molecule using 
TDDFT within the adiabatic approximation. We compared our reliable results with 
experimental data finding a satisfactory agreement even if a part of the excitations 
can not be reproduced by frequency-independent kernels such as the ALDA one. 
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Chapter 7 



A paradigmatic case for 
semicore and spin-polarization 
effects 

in electronic spectra of solids: 
bulk iron 

Semicore states of transition metals such as iron are outside the reach of "standard" 
pseudopotential of DFT which, even when non linear core corrections are adopted, 
include them in the frozen atomic core. In this work we present an analysis of several 
pseudopotentials for iron generated in the Troullier Martins and Hamann scheme 
assuming both local density and generalized gradient approximations (LDA, GGA) 
for the exchange-correlation functional and considering core- valence partitions with 
and without including semicore orbitals among valence states. Non linear core cor- 
rections are considered and the pseudopotential transferability has been checked. We 
calculate structural and electronic properties of the a phase of iron and we present 
a comparison between calculated optical conductivity and experimental data. 

7.1 Pseudopotentials for iron 

Pseudopotentials (PP) are a well-established tool in ab initio structure calculations 
of solids. A review of the topic can be found in literature ranging from the most in- 
fluential works [TU1 HH H21 HH HSl HS1 HU to other important but less fundamental 
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papers [T5J Q2J [20j [21] . Important advantages of the pseudopotential approach can 
be summarized in the following two points: first, by replacing the atom by a pseu- 
doatom, the number of orbitals which have to be calculated is reduced, and, second, 
using plane waves, the size of the basis set can be substantially reduced, because 
the pseudo-wavefunctions are smoother than their all-electron counterparts. 

However in all cases where an overlap between valence and core wave functions 
exists, the frozen core approximation underlying the construction of all pseudopo- 
tentials is not well satisfied. One way to overcome this problem is the inclusion of 
a core correction considering the non linear contribution of the core charge to the 
exchange-correlation potential (NLCC). Another more straightforward solution is 
the explicit inclusion of the semicore electrons into the valence shell. 
In this work we considered both these approaches in order to build pseudopotentials 
for the specific case of iron. Such pseudopotentials are used to calculate the struc- 
tural and electronic properties of the a phase of bulk iron. In conclusion we report 
our results for the optical conductivity compared with experimental data. The 3d 
wavefunctions of iron are strongly localized and show a significant overlap with the 
3s3p orbitals, although the latter are much lower in energy. Looking at Fig. \7.1\ it 
is evident that the separation of the electronic system into well-isolated core and 
valence shells is ambiguous, because of the large overlap between 3s and 3d states. 
For this reason, we considered two electronic configurations: one with 8 (3<i 6 4s 2 
states) and the other with 16 (3s 2 3p 6 3d 6 As 2 states) valence electrons. We used two 
schemes for pseudopotential generation, referring to Hamann [21] and Troullier- 
Martins [20] works respectively. In the former scheme a fixed cutoff radius r c i by 
an exponential function exp[— (r/r c i)A], in the latter scheme wavefunctions are built 
with a parametric form r 1 exp[p(r)] below r c i, where p(r) is a polynomial of order 
six in r 2 . 

In the case of a standard pseudoatom with 8 valence electrons we verified the impor- 
tance of including non linear core corrections (NLCC) to correct GGA pseudopo- 
tentials presenting fake wiggles in the r — > limit due to the gradient dependence 
of the exchange-correlation approximation. 

In the case of a 16 electrons pseudoatom we could not generate a transferable pseu- 
dopotential in generalized gradient approximation, neither in the Troullier-Martins 
nor in the Hamann scheme. The only pseudopotential with good transferability was 
generated within the local density approximation following the Hamann scheme. 
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Figure 7.1: Wave functions of 3s3p3d and 4s level of iron can be distinguished by 
their different localization. The not-negligible overlap between 3d and 3s leads to an 
ambiguity in the definition of core and valence states and to the need of corrections 
to the frozen core approximation. 
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Figure 7.2: Logarithmic derivatives in the case of a 16 electrons pseudoatom. Left: 
Pseudopotentials generated according to the Troullier-Martins scheme show discrep- 
ancies between the differents curves, and are less transferable. Right: Pseudopoten- 
tials generated within the Hamann scheme are more transferables and consequently 
more accurate in reproducing scattering properties. 
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In order to check the pseudopotential transferability we studied logarithmic 
derivatives and compared them to the their all-electrons counterpart. 
Figure 17.21 shows logarithmic derivatives in case of Troullier-Martins (left) and 
Hamann (right) pseudopotentials respectively. Hamann pseudopotential are clearly 
more transferable, and therefore more accurate in describing scattering properties. 
However we should be aware that Hamann pseudopotentials are usually much more 
expensive in terms of plane waves required for convergence than the ones generated 
within the Troullier Martins scheme. As a general rule, depending on the specific 
system and properties we want to describe, a compromise between accuracy and 
efficiency should guide the choice of the suitable pseudopotential. 

7.2 Properties of bulk iron 

In this section we report the results of DFT calculations of structural and electronic 
properties of bulk iron. We also show a Time Dependent DFT calculation of the 
optical conductivity, and we compare our results with experimental data. 

7.2.1 Structural properties 

An useful comparison between the generated PP can be based on the evaluation of 
physical quantities. In particular we calculated the lattice parameter ao of iron. We 
considered iron in his a phase, the most stable at normal temperatures, where the 
metal presents a body centred cubic (BCC) crystal structure. 

The optimized value of ao has been obtained by an estimation of the minimum 
of the Birch-Murnaghan [158] 1159] equation of state for solids: 



ABINIT package [68, 122J to compute converged total energy for different values of 
the cell parameter ao, and we fitted the obtainted points with Eq. 17.11 (see Fig. 17. 
In Tab. l7.ll we report a summary of the lattice parameters and bulk moduli obtained 
with the fitting procedure. Moreover Fig. [73] shows a comparison between LDA and 
GGA results, confirming a well known understimation of the experimental value by 
LDA, that can be partially corrected using GGA. 




(7.1) 



where Bq is the bulk modulus and Vq is the cell volume at the minimum. We used the 
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Figure 7.3: Birch-Murnaghan fit of DFT total energy calculations for iron (8 
valence electrons) performed within local density (continuous line) and generalized 
gradient (dotted line) approximations. Lattice parameter calculated within LDA are 
lower than the experimental value, and GGA partially correct this understimation. 



Pseudopotential 


a 


B 






[Bohr] 


[GPa] 


M 


8e- TM GGA NLCC 


5.56 


67.86 


2.47 


8e- TM LDA NLCC 


5.34 


106.31 


2.11 


16e- H LDA NLCC 


5.27 


101.23 


2.19 


Exp>) 


5.42 


168 


2.22 



Table 7.1: Comparison of lattice parameter for iron calculated with different pseu- 
dopotentials. Exp.( a ) are experimental data reproduced from Ref. [1601 1161] , 

7.3 Electronic properties 
7.3.1 Electronic properties 

We present here our results for the density of states (DOS) around the Fermi en- 
ergy and the electronic band structure including semicore states. These calculations 
have been performed with the ABINIT package [68] in local spin density approx- 
imation (LSDA). We verified that a random sampling of k points in the Brillouin 
Zone (BZ) is more efficient than the use of a regular grid, in fact with 5000 k DOS 
is converged. In Fig. 17.41 two curves represent converged DOS where the two spin 
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channels show an asymmetric distribution allowing to identify the majority and mi- 
nority component. The corresponding band structure is presented in Fig. 17.51 where 
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10 



Figure 7.4: Spin resolved density of states of BCC iron calculated within DFT- 
LDA with ABINIT [BB]. Majority (red line) and minority (blue line) spin channels 
are presented. The Fermi energy is set to zero for clarity. 



high symmetry crystallographic directions are chosen in agreement with previous 
literature [162|, 1163} [L64J. We used Hamann LDA pseudopotentials generated in- 
cluding semicore states in order to describe the core bands 3p and 3s bands. We 
obtain band energies of Fe-3s (85 eV) and Fe-3p (52 eV) as depicted in Fig. 17.51 In 
the same figure we show the effect of the spin on the electronic properties, removing 
most of the accidental degenerancies of band structure and splitting the flat core 
states. In addition we report in Fig. 17.61 (top panel) a comparison between bands 
calculated within LDA or GGA, concluding that the GGA corrections do not change 
dramatically the electronic structure making LDA a satisfactory approximation for 
our scope. Even the use of semicore PP (bottom panel in Fig. 17. 6p do not influence 
considerably the estimation of valence states, therefore it is possible to conclude that 
our approximations are reasonable in order to predict correctly the optical proper- 
ties of bulk iron. 

On the contrary, the case of 3p levels is more subtle. In fact, because the angular 
momentum I is not zero for p levels (i.e. 1=1), the total quantum number j = s + I 
must be considered, and the spin orbit coupling can influence the exact degenerancy 
and position of the bands. Even if this case will not be treated in the present work, 
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state 


A% 


Dispersion 






[eV] 


[meV] 


[eV] 


3p 


3.03 


360 (at N) 




3s 


3.15 


120 (at T) 


4.9 



Table 7.2: Calculated splitting of the core states compared with available experi- 
mental data and calculated dispersion. 

because beyond the aim of the thesis, it is worth to mention that a fully relativistic 
approach is required in order to obtain reliable results for this kind of states. Con- 
versely, the 3s states are correctly predicted within the DFT-LSDA approach and 
the up-down splitting A-m = 3s^ — 3s i is comparable with the available experimen- 
tal data [161J. For completeness, we report in Tab. 17.21 the calculated splitting and 
dispersion of both the 3s and 3p levels. 

In conclusion we also report an analysis of the relation between total magnetiza- 
tion \i of iron and the 3s splitting. It was suggested in the past j!65[ 1166] that s core- 
level splitting could be used to monitor the magnetic moment or the hyperfine field, 
because the splitting should vary linearly with the spin state of the unfilled inner 
shell. X-ray photoelectron spectroscopy (XPS) measurements for nonmetallic tran- 
sition metal compounds, on rare-earth metals and ionic compounds [165, 11661 H67] 
are compatible with this scheme. 

However, Fe-3s XPS splitting studied in crystalline and amorphous alloys have 
shown a poor correlation between the 3s splitting and the magnetic moment of 
the solid (see Fig. 3 in Ref. |161] ) . Our calculations show a linear relation of such a 
splitting with respect to the total magnetization of bulk iron (see Fig. l7.7l top panel), 
giving some basis to better investigate the apparent disagreement with experimental 
data. 

We calculated the total magnetization defined as the difference between the ma- 
jority and minority spin density integrated over the unit cell /i = J n P^(r) — Pi(r). 
This quantities has been evaluated for several lattice parameter values of bulk Fe, 
centred around the experimental one (5.42 Bohr). Changes of the lattice parameter 
influence the local distance between atoms, reproducing the effect of a change of 
pressure. Moreover changes on lattice geometry affect total magnetization. The 
dependence of the splitting with respect to the pressure is shown in Fig. 17.71 bottom 
panel. 
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Figure 7.5: Spin resolved band structure of BCC iron along standard high sim- 
metry directions of the Brillouin zone, calculated in LDA approximation with the 
ABINIT [68J code. Valence (top) and semicore (bottom) states are distinguished in 
different panels in order to evidence energy scales. 
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Figure 7.6: Spin resolved band structure (valence states) of BCC iron calculated at 
the experimental lattice paramter. Top: comparaison between GGA and LDA calcu- 
lation; Bottom: comparaison between LDA calculations with and without semicore 
states. 
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Figure 7.7: Splitting of 3s levels as a function of the lattice parameter (top panel) 
and the magnetization (bottom panel). 
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7.4 Optical conductivity 

In this last section we report the calculations of the dielectric function (e = e\ + 
%£•}) obtained within the indipendent particle random phase approximation (IP- 
RPA) and including crystal local fields effects (RPA-LF). This quantity is closely 
related to the electronic structure and follows straightforward from theory, however 
experiments are often presented in terms of the optical conductivity (a = u\ + io"2, 
see Refs. [168] ). However, e and a are connected by the well known relations: 

- = £ ™ 

In Fig. 17.81 we show the optical conductivity of bulk iron calculated with the DP 
code [78] . There is a nice agreement between the calculated curves and the experi- 
mental data reported in Ref. [168J with respect to the general shape and the position 
of the maxima. We have not included any intraband contribution in our calculation 
and this is the reason why the sharp structure in the low-energy region of the ex- 
perimental spectra is not reproduced (Fig. 17. 8p . We report, for comparison, a recent 
calculation of the optical conductivity including the Drude peak (see Ref. JT69J for 
details about the treatment of the intraband transitions), in order to show that only 
the low energy part of the spectum is affected by this approximation. 
Spin flip is also not allowed and this bring us to conclude that the maximum of 
the computed conductivity near 2.5 eV results from transitions between collinear 
spin states, in particular between states of the majority spin channel. It is easy to 
identify these excitations looking at Fig. 17. 4t transitions connect occupied states 
and the unoccupied states just above the Fermi energy for the same spin channel 
(dashed line). Crystal local fields have negligible effects, and this is in agreement 
with the cases of others metals. As shown in Fig. 17.81 bv the red and blue lines. 
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Energy [eV] 

Figure 7.8: Comparison between optical conductivity experimental data from 
Refs. |168| and RPA calculation performed with the DP code [78]. Results with 
and without the inclusion of crystal local field effects (LF) are shown by the red and 
blue lines, respectively. 

7.5 Conclusions 

We succeeded in generating a set of pseudopotentials for iron, a prototypical tran- 
sition metal for which an obvious separation between core and valence orbitals does 
not exist. 

Semicore states have been either explicitly included in the valence set, or treated 
within the usual non linear core corrections scheme. Our set of PP has been tested 
against structural and electronic properties of BCC bulk iron, including the lattice 
constant, spin resolved density of states, band structure. 

We also reported an analysis of the Fe-3s splitting with respect to the total 
magnetization. We found a linear dependence between these two quantities and 
between splitting and the lattice parameter. This results support the use of the 
s core-level splitting as a monitor of magnetization, and are compatible with a 
similar behaviour presented by non metallic transition metal compounds, rare earth 
metals and ionic compounds. Moreover we found that the 3s splitting calculated at 
the experimental value of the lattice parameter substantially agrees with XPS core 



7.5. CONCLUSIONS 



127 



level measurements reported in Ref. [161] . However the linear trend theoretically 
predicted is in apparent disagreement with experimental results where 3s XPS peak 
splittings do not correlate with the Fe magnetic moments. 

Concerning the optical conductivity good agreement between RPA results and 
the experimental data is found. This allows one to conclude that the main spectral 
feature are well reproduced by spin-collinear excitations. 
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Chapter 8 

Iron, cobalt and nickel pyrites 



Transition-metal pyrites form a series of compounds with a large variety in elec- 
trical, magnetic and optical properties. In particular, recently there has been a 
renewed interest in iron, cobalt and nickel disulfides because of their potential in 
future technological applications. Iron and nickel-controlled doping of C0S2 gives 
rise to a tunable source of highly spin-polarized electrons |170j . This property is 
extremely interesting to design new devices exploiting the spin character of the elec- 
trons in addition to their charge. In fact the essence of the current focus area termed 
spintronic, or spin- electronics, is to use the electron's spin, as well as its charge in 
creating new devices or enhancing the functionality of the existing ones [171} 1172] . 

In this section we will present structural, magnetic and electronic properties of 
FeS2, C0S2 and N1S2 calculated within the density functional theory framework. A 
comparison of our results with respect to experimental data and to previous calcu- 
lations is also discussed. 

8.1 Motivations 

The rapid developement of spin valve-based magnetic read heads and the emergence 
of spintronic [171} 1172] has thrown up a need for a better understanding of spin- 
polarized materials [173J. 

Spintronic is based on the up or down spin of the carriers rather than on electrons 
or holes as in traditional semiconductors electronics. In particular spin-polarized 
transport will occur naturally in any material for which there is an imbalance of 
the spin populations at the Fermi level. This imbalance is present in ferromagnetic 
metals where the density of states for spin up and down electrons are shifted in 
energy with respect to each other. Commonly in these materials there is an unequal 
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Figure 8.1: Schematic representation of the density of states in case of a metal 
(left) and an half metal (right). 

filling of the bands, which is the source of the net magnetic moment, and causing 
the up and down carriers at the Fermi level to be different in number, character, 
and mobility. This inequality produces a net spin polarization but the sign and 
magnitude of that polarization depends on the specific material. Materials having 
at the Fermi level only one occupied spin band are usually called half metals (see 
Fig. EH). 

A fundamental component in any spintronic device is a ferromagnetic electrode 
which is used as a source of polarized electrons, and an high value of spin polariza- 
tion P=(N^-N|)/(N|+N|) at the Fermi level can provide significant benefits. When 
the spin polarized electron current crosses a sample having a non zero total magneti- 
zation the only states that are available to the carriers are those for which the spins 
of the carriers are parallel to the spin direction of those states at the Fermi level. If 
the magnetization of the material is reversed, the spin direction of those states also 
reverses. Thus, depending on the direction of magnetization of a material relative 
to the spin polarization of the current, a material can be either a conductor or an 
insulator for electrons of a specific spin polarization (see Fig. 18. 2p . 

The largest effect is generally seen for the most highly polarized currents, there- 
fore, there are continuing efforts to find 100% spin-polarized conducting materi- 
als. However, partially polarized materials (such as Fe, Co, Ni and their alloys), 
are adequate to develop technologically useful devices and can show unexpected 
properties. For example, transition metal compounds have beed attracting ex- 
tensive attention for that reason. Among them, FeS2, C0S2 and N1S2 have been 
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Figure 8.2: Schematic rep- 
resentation of spin polarized 
transport from a FM material 
to a metal and a FM mate- 
rial again. According to the 
reciprocal configuration of the 
spins, transport is allowed or 
forbidden. 



studied from the experimental pH UM Eg [T7T\ and theoretical [ESI LTZSl ESQ] 
point of view. Moreover these compounds have been used to build and to study 
Coi_ a .Fe :r S2 and Nii_ x Co a; S2 compounds, recently predicted to be tunable half met- 
als (see Refs. HSU 1771 [HH [L83l MM for a review of works on Coi__ a .Fe x S 2 and 
Ref. [ISH] for Ni^Co^Ss). 

In the following sections we will focus on structural and magnetic properties of 
FeS2 and C0S2, as a first step to be able in the future to describe more complex 
compounds. 

8.2 Structural and magnetic properties 

All transition metal disulfides (MS2 with M a 3d-transition metal atom), crystallize 
in a cubic pyrite structure of space group T®(Pa3) in which metal atoms are located 
in face-centred postions. Structure can be considered as an NaCl-like grouping of 
metal and chalcogen atom pairs (sulfurs). In Fig. 18.31 the atoms arrangement is 
shown from three perpendicular views. Every metal atoms (dark grey circles) is 
surrounded by six nearest-neighbour sulfurs in a distorted octahedral environment, 
while each sulfur (light grey circles) bonds to one sulfur (S-S bond) and three metals 
in a distorted tetrahedral arrangement. Distance of sulfur-sulfur pairs (S-S) is short 
because of a covalent bond. The formation of S-S pairs is characteristic feature of 
these structures. 

In particular metal atoms are located at the positions (0,0,0), (0,1/2,1/2), (1/2,0,1/2) 
and (1/2,1/2,0), the eight sulfur atoms instead are located at position ±(u,u,u), 
±(u + 1/2,1/2 -u,u), ±(u,u + l/2, 1/2 -u), ±(1/2 - u,u,u + 1/2). The values of 
u and a (the lattice parameter) are taken from Wyckoff [186J , in particular we used 
u = 0.386 and a = 5.407 A in case of FeS2, u = 0.389 and a = 5.524 A in case of 
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Figure 8.3: Atoms rearrangement in pyrite structure along three orthogonal di- 
rections. Dark grey circles represent metal atoms, smaller light grey circles depict 
sulfur atoms. Graphs are generated with Xcrysden package |136j . 



System 


a [A] 
Ref. [186] 


s-s [A] 

this work 


s-s [A] 

from Ref. [l78j 


fjt/at. [hb] 
this work 


///at. [fi B ] 
Ref. [17411181] 


FeS 2 


5.407 


2.18 


2.14 


~ 


~ 


CoS 2 


5.524 


2.29 


2.12 


0.98 


0.9 


NiS 2 


5.677 


2.06 


2.06 








Table 8.1: Optimized S-S distances calculated for FeS 2 , CoS 2 and NiS 2 are com- 
pared with values reported in Ref. [178j . The calculated ground state magnetization 
per metal atom is compared to experimental values reproduced from Ref. [174"[ 118 J] 

CoS 2 and u = 0.395 and a = 5.677 A in case of NiS 2 . 

We performed a geometric optimization using the BFGS algorithm [141|, 11421 
H331 ESI with the ABINIT [68] code. The optimized S-S distances of FeS 2 , CoS 2 and 
NiS 2 are reported in Tab. 18. 1\ and compared with values reported in Ref. [178j . All 
the ground state calculations have been performed within the DFT-GGA frame- 
work (Perdew-Burke-Ernzerhof (PBE) parametrization with a cutoff energy 
E cu t = 34 Ha. The geometry relaxation is performed by setting a tolerance of 
0.02 on the ratio of differences of forces to maximum force, reached twice succes- 
sively, will cause a self consistent cycle to stop. The agreement of the computed S-S 
distance with respect to previous calculations is reasonable, as listed in Tab. 18.11 
Moreover, Figure 15741 shows the electronic density distribution for a suitable value of 
the chosen isosurfaces in order to evidence the formation of the characteristic S-S 
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Figure 8.4: Iso- 
surfaces of electronic 
density (grey regions) 
around the atoms of 
CoS2- Cobalt atoms 
(light blue) and sul- 
fur atoms (yellow) are 
represented inside the 
crystal structure. The 
point of view is cho- 
sen in order to evi- 
dence the characteristic 
S-S bond in the center 
of the crystal (red ar- 
rows). Graphs are gen- 
erated with the Xcrys- 
den package [136J. 

bond in the center of the crystal structure. 

The magnetic character of the electronic ground state of all the compounds has 
been also analyzed. The values of total magnetization per metal atom are reported 
in Tab. 18. 11 with experimental reference values reproduced from literature. 
The three compounds display different behaviours: FeS2 presents a total magnetiza- 
tion close to zero, C0S2 has a ferromagnetic ground state with total magnetization 
per cobalt atom close to 1 /xb, and N1S2 is non-magnetic. Experimental studies on 
transition metal disulfides (see Ref. [174] ) revealed the non magnetic nature of all 
this class of compounds, except for the case of CoS2- Experiments indicate that FeS2 
as well as N1S2 are paramagnetic semiconductors. In addition, at low temperature 
N1S2 presents a transition to an antiferromagnetic phase. 

Our calculations predict correctly this behaviour. In particular for C0S2 we find a 
ground state electronic configuration with a total magnetic moment /j, =0.98 per 
cobalt atom at T=0° K, which is close to the experimental value /i =0.9 \ib (see 
Ref. [Trap . 

In addition we report in Fig. I8.5l the electronic density isosurfaces that evidences 
the distribution of spin up and spin down components. In the case of FeS2 the spin 
up and down charge distributions are located at the same positions without any 
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motif distinguishig between the two components. On the contrary, in case of C0S2 
it is possible to appreciate a regular ordering of the electronic density distribution 
according to the two spin components (red and blue respectively in Fig. 18. 5 [) . The 
maxima of the distributions are located at the metal sites of the crystal and the 
up-down spin density is alternated in adiacent sites. 

NIS2 still presents a regular ordering of the spin density distributions but, as the 
total magnetization is found to be zero, suggests an antiferromagnetic character. 
Even if the study of MS2 is beyond the scope of this thesis, we notice that this 
result is in agreement with Rcf. [174J. 

FeSL Col NiS, 




Figure 8.5: Electronic density isosurfaces for the up (red) and down (blue) spin 
component for FeS2 (left), C0S2 (center) and MS2 (right). Graphs are generated 
with the Xcrysden package [136J. 



8.3 Electronic properties 

In order to analyze the electronic properties of FeS2 , C0S2 and N1S2 we calculated 
the spin resolved density of state (DOS) within DFT-GGA (PBE parametrization). 
The software ABINIT [68] has been used in order to perform convergence study 
and calculate the final distribution. For F3S2 we used a Monkorst Pack [83] grid of 
12 x 12 x 12 k points in the Brillouin zone (BZ), corresponding to 76 k, applying 
a rigid shift of 0.5 in the three spatial directions. For C0S2 the grid consisted of 
10 x 10 x 10 k points and 6x6x6 for NiS2- A thermal broadening (cold smearing) 
tsmear = 0.007 Ha is applied to simulate the metallic occupation of levels following 
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the recipe of Ref. [187^1 . 

Density of states are calculated using the tetrahedron method giving a faster conver- 
gence over the number of k poins required with respect to the standard integration 
over the Brillouin zone. 

In Fig. 18.61 the DOS of FeS2, C0S2 and MS2 are shown. The three compounds 
present a similar density of states but different electronic behaviour. We discuss now 
our results that are in general agreement with experimental works [1751 11761 1188] 
and theoretical calculations based on linear combination of atomic orbital (LCAO) 
or semiempirical and self-consistent tight-binding (TB) approach |178} 1179} 1189] . 
The first two lowest energy bands are associated with bonding and antibonding pairs 
of orbitals of the sulfur atom dimers (sa and scr* orbitals). The following complex 
structure is due to sulfur 3p and the metal 3d orbitals. Then the crystal field of 
the disulfide anions (sulfurs) splits the cation (metal) 3d non-bonding orbitals into 
three low-lying 3d(t2j) and two higher-energy 3d(e g ) levels (see Fig. 18. 6|) . 
All the disulfides studied have completely filled 3d(t2j) orbitals, and as the atomic 
number increases the additional electrons (none for iron, one for cobalt and two for 
nickel) fill in the 3d{e g ) orbitals. The different filling of these bands determine the 
electronic nature of the three compounds studied. 

In particular, our calculations predict FeS2 as a small band gap semiconductor, 
in agreement with experiments (optical and conductivity measurements [176 J and 
X-ray photoemission spectoscopy [175} 1190] ) and theoretical works [1781 1179j . From 
our calculations the value of the FeS2 band gap turns out to be 1.2 eV and 0.98 eV 
for the spin-up and spin-down components respectively (Fig. 18.61 top panel). Even 
if a more detailed analysis including band structure calculations is required, we can 
already conclude that the calculated values are close to the experimental ones rang- 
ing from 0.9 eV to 1.2 eV (see Ref. |176| ) and for which there is no general consensus. 
Difficulties in the exact experimental determination of the band gap of FeS2 could 
rise because of a large difference between the indirect and direct gap as pointed out 
by Zhao et al. |179j . In fact in that work all calculations in the local density ap- 
proximation lead to a minimum indirect gap Ej 9 = 0.59 eV and to a smallest direct 
gap E dg = 0.74 eV. 

Concerning C0S2, we found similar structures for low-energy states with respect 
to FeS2, but shifted to lower energies (see Fig. 18. 6\ central panel). Moreover our 
DFT-GGA calculations confirm the half metallic nature of C0S2, where the two spin 



lr The smeared delta function is denned by: 5 S = -j^ e x ( axZ ~ — |aa; + |) where a = 
—0.5634, this choice minimizes the bump. 
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components present density of states characteristic of a semiconductor (red line) or 
a metal (blue line) . In Fig. 18.61 the Fermi energy is set to zero for clarity. 
It is worth to mention that half-metallicity of C0S2 has been recently discussed in 
literature from both the experimental and theoretical points of view. Direct mea- 
surements of spin polarization and magnetotransport [170] show that C0S2 is not 
completely polarized (i.e. P = ^(g^ j-ivjxEf) * s no ^ ^ ar S e )- On the contrary, the 
temperature variation of the electronic structure, studied by means of optical re- 
flectivity measurements confirms the half metallicity of the compound. From the 
theoretical point of view Shishidou et al. [180J presented a detailed comparison be- 
tween LDA and GGA electronic structure calculations within density functional, full 
potential linearized augmented plane waves calculations. In that work the authors 
conclude that GGA greatly modifies the LDA band structure from metallic to half 
metallic. These conclusions are in agreement with our finding displayed in Fig. 18.61 
In the case of MS2 the spin density of states presents a symmetric behaviour that 
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Figure 8.6: Spin resolved density of states of FeS2, C0S2, MS2 calculated with 
tethraedron method with the software ABINIT [68]. 



confirm the apparently non magnetic nature of that compound and the zero total 



8.4. CONCLUSIONS 



137 



magnetization obtained from our calculations (see sec, \8.2\i . However, as underlined 
preivously, N1S2 presents a local regular ordering of the spin density distribution jus- 
tifying its antiferromagnetic nature (see Refs. |191j ). Among the transition metal 
disulfides, N1S2 is considered to be a Mott-Hubbard insulator [192J. For that reason 
more specific theories are required in order to correctly describe N1S2 properties. 

8.4 Conclusions 

We presented a systematic study of structural, magnetic and electronic properties 
of FeS2, C0S2 and N1S2, three examples of 3d-transition metal sulphides with cubic 
pyrite structure. As anticipated in the introduction of this chapter, the interest on 
these systems is recenlty renewed, due to the great potential of such systems for 
spintronic applications. 

From our calculations we obtained the S-S relaxed distance, characteristic of 
these compounds, the ground state total magnetization and the spin resolved den- 
sity of states. We correctly predict that C0S2 is the sole magnetic system of this 
class of compounds, with a total magnetization per metal atom \x =0.98//b, in agree- 
ment with experimental values. On the contrary FeS2 and MS2 do not show any 
total magnetization, the former being non-magnetic and the latter presenting lo- 
cally a regular arrangement of the two spin components of the electronic density 
distribution, typical of an antiferromagnetic material. Moreover DOS calculated for 
C0S2 is typical of an half metal, while FeS2 and NiS2 typical of small band gap 
semiconductors . 

All these results are in nice agreement with experimental data and previous theo- 
retical calculations. 

These results now open the way to study also the optical properties of such com- 
pounds. 

In conclusion, as anticipated in the introduction of this chapter, the interest on 
these systems is recenlty renewed due to the great potential for the emerging spin- 
tronic technological applications. In fact, when C0S2 is doped with iron or nickel 
in a solid solution like Coi_ a; Fe a: S2 or Coi_ x Ni :r S2. the metal to semiconduc- 
tor transition can be tuned increasing the iron or nichel concentration. More- 
over the spin polarization can be controlled generating an highly polarized electron 
source [1771 HHH M [[831 H831H55] . 
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Conclusions 



This thesis is devoted to ab initio calculations of ground state and excited state 
properties of realistic systems within the density functional theory (DFT) and its 
Time Dependent generalisation (TDDFT). 

We used theoretical spectroscopy tools in order to study several systems with 
different dimensionality (surfaces, molecules, bulk crystals). Starting with the di- 
electric function e(cu), obtained by the response function in linear regime, we cal- 
culated the anisotropy reflectivity (RA) spectra and the reflectance electron energy 
loss (REEL) spectra for the Si(100) clean and oxidized surfaces. In the case of the 
clean surface, we considered three surface reconstructions p(2 x 1), p{2 x 2) and 
c(4 x 2). Thanks to the comparison between experiments and numerical simulation, 
we were able to rule out the p(2 x 1) reconstruction and to define the origin of the 
REEL peaks without ambiguity. 

In the second part of the work, we evidenced the problem of the correct description 
of excitation spectra of open shell molecules within DFT-LDA. We calculated the 
energy levels and the first excitation energies for the BeH molecule in the TDDFT 
framework. These calculations have been a first step in order to approach the case 
of magnetic metals and half metals. In fact, this part of the work was dedicated 
to the study of optical properties of magnetic iron alloys, that are very interesting 
materials for new applications in the spintronic domain. 

In this framework we evaluated ground state properties and conductivity of the bec 
crystal phase of iron in order to validate the theoretical approach comparing the 
results with experimental data. 

Finally we studied the ground state properties of some 3<i-transition metal disul- 
fides, i.e. FeS2, C0S2 and MS2, having a commun cubic pyrite structure but different 
electronic and magnetic properties. Starting from these compounds it is possible 
to create more complex systems such as iron or nickel doped C0S2 alloys, that are 
interesting for the design of new technological devices based on spin electronics. 
From the numerical point of view, we implemented an original method in the 
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large scale ab initio code DP in order to calculate the indipendent particle dynami- 
cal response function x°( r > r '> w )> built from the eigenvalues and eigenvectors of the 
Kohn and Sham hamiltonian. We have demostrated that the method, based on 
the Hilbert transform, is efficient for large size systems, as in the case of surfaces. 
Moreover, we have generalized the code to spin variable in order to study magnetic 
properties of realistic applications. 



Milano, December 2008 
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Appendix A 



Determination of second 



derivative spectra 



We implemented a Savitski-Golay (SG) smoothing algorithm [193J in order to com- 
pute the derivative spectra. 

The SG algorithm fits the EEL curve with polynomials preserving spectral features: 



n=-n L 

where gi represents the smoothed function, and cn the coefficients obtained by the 
coefficient matrix of polynomial of degree M by the relation: 



where Aij = V is the coefficient matrix of polynomial of degree M and e n is the 
unit vector. Then SG algorithm returns the derivative of fc-degree of the smoothed 
curve. When k > 1 the coefficients in equation lA.il have to be normalized replacing 
cat with cjvfc!. 

Smoothing step is required to correctly identify the physical EEL peaks, in fact, 
derivative can dramatically enhance the noise of the EEL spectra due to the discrete 
k point sampling. See Fig. lA.ll for an example of the method. 




(A.l) 



cat = (A T A)-\A T e n ) 



(A.2) 
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Eriuiyy (t^VJ 

Figure A.l: Computation of the second derivative spectra (blue line) of a typical 
EEL spectrum. The bare EEL spectrum (black line) and the resulting fit with 
polynomials (red line) are also shown. 
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